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CHAPTER  1 
INTRODUCTION 

Radar,  or  itodio  /Electing  and  hanging,  has  been  used  for  decades  to  detect  distant  objects, 
measure  their  velocity,  and  more  recently,  for  terrain  mapping.  It  was  theorized  in  the  early 
1900s.  and  first  demonstrated  in  the  mid-30s.  Subsequently,  it  was  developed  to  detect  airplane 
bombers  during  WWII  [l] .  Since  then,  it  has  been  refined  and  used  in  hundreds  of  applications 
such  as  space  exploration,  airplane  avoidance  systems,  insect  and  bird  tracking,  weather  obser¬ 
vation  and  prediction,  and.  of  course,  highway  speed  enforcement. 

When  radar  is  used  to  simply  detect  the  range  to  a  distant  object,  a  series  of  electromag¬ 
netic  pulses  is  transmitted  in  the  direction  of  the  object  and  the  two-way  time  delay  or  the 
reflected  pulses  is  measured.  Knowing  the  pulses  propagate  at  the  speed  of  light,  we  can  readily 
calculate  the  range  to  the  object.  If  we  wish  to  measure  the  radial  velocity  of  an  object,  we  can 
transmit  a  continuous  wave  (CW)  electromagnetic  signal  (a  simple  sinusoid)  and  measure  the 
shift  in  frequency  of  the  returned  signal  induced  by  the  moving  object.  Again,  it  is  simple  to 
calculate  the  velocity  of  the  object  from  the  frequency  shift  (Doppler  shift)  and  the  known 
propagation  velocity  of  the  signal. 

When  the  azimuthal  position  (perpendicular  to  range)  is  also  desired,  the  transmitted  sig¬ 
nal  must  be  sent  as  a  narrow  beam  which  is  no  wider  than  the  desired  azimuthal  resolution.  It 
is  easily  shown  that  at  a  range  R  and  carrier  frequency  wc,  an  antenna  with  an  aperture  size  D 
will  have  an  azimuthal  resolution  of  order  Ro>c/D  meters.  For  radar  operating  at  conventional 
microwave  frequencies  (10’  Hz),  this  implies  that  very  small  resolution  requirements  dictate  a 
physical  antenna  of  gigantic  proportions  -  several  thousand  feel  across!  For  an  airborne 
antenna,  this  is  impractical;  however,  very  high  azimuthal  resolution  can  be  achieved  by  syn¬ 
thesizing  a  large  phased  antenna  array  in  an  imaging  system  called  synthetic  aperture  radar. 

Synthetic  aperture  radar  (SAR)  is  a  means  of  obtaining  high  resolution  terrain  maps  (com¬ 
plex  reflectivity  maps)  by  coherently  processing  the  backscattered  radar  signal’s  phase  (referred 
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to  as  “phase  history  ’).  In  1953.  the  University  of  Illinois  experimentally  demonstrated  strip 
mapped  SAR[2}.  This  was  the  earliest  form  of  SAR  in  which  the  airborne  radar  antenna  is  held 
at  a  fixed  squint  angle  relative  to  the  terrain  patch  of  interest.  Coherent  processing  of  the 
returned  radar  signal  produces  an  effective  aperture  many  times  the  size  of  the  physical 
antenna,  thus  creating  an  azimuthal  beamwidth  proportionally  narrower,  resulting  in  finer 
azimuthal  resolution. 

A  higher  resolution  form  of  SAR  is  called  spotlight  mode  SAR.  This  mode  requires  that 
the  antenna  be  continuously  positioned,  or  steered,  toward  the  center  of  the  terrain  patch  dur¬ 
ing  the  plane's  flight.  The  same  patch  is  effectively  pictured  from  many  different  angles,  and 
when  coherently  processed,  produces  a  much  higher  resolution  image  than  could  be  obtained 
from  a  single  observation  angle.. 

Initially.  SAR  reflectivity  data  were  processed  optically.  The  phase  histories  were  recorded 
on  film1  and  then  processed  on  an  optical  bench  through  a  complex  and  bulky  arrangement  of 
lenses  and  coherent  light  [3. 4] .  The  data  w ere  continuous  in  the  range  direction  and  sampled 
in  the  azimuth  coordinate  (the  discrete  return  pulses).  An  optical  Fourier  transform  is  thus 
continuous  in  range  and  discrete  in  azimuth  and  computed  almost  instantaneously.  Although 
high  resolution  was  achieved,  the  characteristics  of  the  optical  processing  equipment  made  the 
reconstruction  procedure  cumbersome  in  the  laboratory  and  impossible  during  inflight  data  col¬ 
lection.  The  advent  of  very  fast  digital  processors  and  inexpensive  semiconductor  memory 
made  real-time  imaging  feasible,  i.e..  the  terrain  maps  could  be  (theoretically)  generated  in  the 
data  gathering  platform  online  and  displayed  during  the  flight.  Digital  processing  of  the  signals 
led  to  other  problems,  however,  and  such  issues  as  polar  sampling  criterion  and  polar-lo-digital 
interpolation  became  the  focus  of  much  study.  Stark  [5]  presents  some  polar  sampling 
theorems  for  the  case  of  complete  2ir  sampling  rasters  such  as  those  used  in  computer-aided 
tomography  (CAT):  however,  they  are  difficult  to  apply  to  the  case  of  spotlight  mode  SAR 

'A  light  source,  modulated  by  the  return  signal,  is  mechanically  moved  over  tilm  at  a  speed  proportional  to  the 
sweep  rate  of  the  chirp  radar. 
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where  the  recorded  data  region  is  a  small  slice  of  a  toroid  instead  of  the  complete  disk  shaped 
region  of  CAT. 

In  the  polar  format  algorithm  of  recording  the  phase  histories  of  the  SAR  data,  the  sam¬ 
pled  data  (after  coherent  demodulation)  are  placed  on  a  small  section  of  a  polar  grid.  If  an 

inverse  discrete  Fourier  transform  (IDFT)  for  polar  grids  were  available,  then  the  radar  data 

i 

could  theoretically  be  inverse-transformed  directly  and  the  final  image  displayed  on  a  polar 
|  display  device.  Such  an  IDFT  is  not  known,  however,  and  the  data  must  be  interpolated  to  a 

rectangular  grid  before  the  IDFT  operation.  Hven  if  such  an  IDFT  were  possible,  polar  displays 
I  are  not  readily  available  (  though  work  is  being  done  in  this  area  [6] ).  and  so  interpolation  in 

the  spatial  domain  would  be  required  for  display  on  a  rectangular  raster.  There  already  exist 
several  fast  Fourier  transform  (FFT)  algorithms  as  well  as  special  purpose  FFT  hardware  avail¬ 
able  to  perform  the  IDFT  for  rectangular  sampled  data  arrays. 

The  most  computationally  intense  task  in  generating  the  radar  image  involves  two- 

i 

dimensional  interpolation  from  the  polar  grid  to  the  rectangular  grid.  While  an  N  by  N  2-D 
FFT  requires  0(N2log2N)  operations,  the  2-D  interpolation  step  often  requires  <XK*N2)  where 
K  is  a  number  ranging  from  logjN  to  as  high  as  N4.  depending  on  the  interpolation  algorithm 
and  order.  Thus,  the  inur^olation  step  enormously  overshadows  the  FFT  step  in  computa¬ 
tional  complexity,  i.e..  processing  time.  This  makes  the  SAR  imaging  tool  difficult  to  implement 
in  real  time. 

The  interpolation  requirements  of  SAR  are  very  similar  to  those  of  the  direct  Fourier 
reconstruction  techniques  used  in  CAT.  Jenkins,  Munson  and  0'Brien(7j  have  developed  an 
analogy  between  SAR  and  CAT  which  places  projection  data  on  the  polar  grid.  There  are  two 
main  differences  between  SAR  and  CAT: 

(1)  SAR  uses  coherent  imaging  and  retains  the  phase  information  from  its  complex  targets 
while  CAT  is  non-coherent  (this  is  not  true  of  diffraction  tomography  which  is  a  coherent 
>  system  that  records  magnitude  and  phase)  and  thus  only  the  magnitude  of  the  projection 
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is  recorded.  The  recorded  CAT  data  are  also  strictly  positive,  while  SAR  data  are  complex 
and  may  have  negative  real  and  imaginary  parts. 

(2)  CAT  Fourier  data  lie  in  a  polar  grid  which  usually  spans  the  full  2 it  in  angle,  and  radially 
from  0  to  some  rm„.  while  SAR  data  lie  in  a  polar  annulus  whose  width  is  proportional  to 
the  radar  signal  bandwidth,  and  which  is  far  removed  from  the  origin  occupying  only  a 
few  degrees  of  angular  width  (typically  2  to  10  degrees). 

For  typical  systems,  the  geometry  of  the  SAR  data  grid  is  very  close  to  a  rectangular,  or  a 
rotated  rectangular  grid.  The  polar  grid  in  CAT  does  not  share  this  property.  Given  the  small 
data  record  of  collected  phase  histories,  and  the  fact  that  the  recorded  Fourier  section  is  oiTset 
from  the  Fourier  origin  (no  dc  in  azimuth),  it  is  surprising  that  one  can  obtain  any  azimuthal 
resolution  at  all.  Munson  and  Sanz  [8]  have  demonstrated  that  it  is  the  coherency  of  the  SAR 
system,  together  with  the  complex  nature  of  the  "spatial"  targets,  which  allows  such  high 
image  quality  from  such  a  limited  amount  of  Fourier  data.  Despite  these  major  differences 
between  SAR  and  CAT  imaging,  many  of  the  digital  processing  techniques  from  CAT  can  be 
applied  to  SAR.  Modifications  to  these  algorithms  are  made  in  accordance  to  the  geometric 
difference  between  them  and  some  approximations  can  be  made  in  SAR  which  are  inappropriate 
for  tomography,  but  the  basic  results  of  sampling  and  windowing,  as  well  as  reconstruction, 
and  error  analysis  of  the  computed  tomography  problem  are  applicable  to  SAR. 

The  reason  that  the  interpolation  requirements  are  so  severe  is  due  to  the  nature  of  the 
transform  operation.  Every  point  in  the  Fourier  domain  contributes  to  the  spatial  reconstruc¬ 
tion.  Single  errors  in  the  Fourier  domain,  therefore  affect  the  entire  spatial  domain.  It  is  very 
important,  then,  for  the  interpolator  to  be  as  accurate  as  possible,  especially  in  the  Fourier  areas 
where  most  of  ihe  energy  is  concentrated.  Without  a  priori  information  about  the  spectral 
energy  distribution,  it  is  thus  imperative  to  be  uniformly  accurate  throughout  the  region  We 
shall  see  that  there  is  a  trade-off  between  interpolator  accuracy  (and  thus  image  resolution)  and 
computation.  An  important  part  of  this  research  is  dedicated  to  reducing  the  latter  two  while 
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maintaining  the  former. 

Multidimensional  interpolation  is  also  important  in  other  problems  where  the  data  gather¬ 
ing  methods  provide  non-rectangular  sampling  geometries,  e.g..  underwater  mapping  [9], 
seismic  mapping  [9],  imaging  the  ocean  and  sea-ice  {10}.  and  three-dimensional  medical  imag¬ 
ing  111.  12. 13]  all  gather  data  on  sampling  grids  which  may  not  be  rectangular.  Digital  rota¬ 
tion.  magnification,  and  reduction  of  images  may  also  require  two-dimensional  interpolation. 

The  locus  of  this  thesis  is  on  the  effects  that  various  two  dimensional  interpolation 
schemes  in  the  Fourier  domain  SAR  data  have  on  spatial  domain  reconstructions.  A  review  of 
the  fundamental  spotlight  mode  SAR  equations  is  presented  in  Chapter  2.  The  tomographic 
formulation  of  SAR  is  also  reviewed  and  compared  to  the  radar  Doppler  interpretation.  Alter¬ 
nate  inversion  techniques  are  also  briefly  mentioned.  A  more  detailed  discussion  of  the  interpo¬ 
lation  problem  in  given  in  Chapter  3.  Classical  techniques  are  reviewed  and  extended  to  the 
SAR  geometry.  Spline  interpolants  are  then  presented  from  a  DSP  point  of  view,  and  a 
clarification  of  the  term  splint  is  given.  The  topics  of  windows,  aliasing,  spurious  targets  and 
artifacts  are  discussed  in  Chapter  4.  as  well  the  formulation  of  alternate  sampling  rasters.  A 
rather  lengthy  set  of  empirical  results  for  many  different  interpolator  algorithms  and  sampling 
strategies  is  given  in  Chapter  5.  It  was  felt  that  more  insight  could  be  obtained  by  studying 
these  reconstructions  rather  than  by  producing  complicated  expressions  for  the  spatially  vary¬ 
ing  polar  interpolators.  Finally,  a  summary  of  the  ideas  and  results  presented  in  this  work  is 
presented  in  Chapter  6.  Topics  of  further  research  to  be  done  in  the  area  of  digital  SAR  image 
reconstruction  is  also  discussed. 
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CHAPTER  2 

SYNTHETIC  APERTURE  RADAR:  THEORY  AND  BACKGROUND 

In  order  to  understand  the  interpolation  problem,  it  is  useful  to  see  how  the  specific 
geometry  developed.  In  the  following  sections,  the  SAR  equations  will  be  reviewed  and  the 
polar/rectangular  geometry  will  be  studied.  The  relationship  to  tomography  will  also  be  exam¬ 
ined. 

2.1  Derivation  of  Exact  Point  Target  Response 

Schwartz  [14]  investigated  three  data  formats  for  SAR  and  included  much  of  the  analysis 
that  was  developed  from  Weis  and  Jenkins  at  Lockheed  [15].  The  analysis  that  follows  was 
the  basis  for  the  original  SAR  simulation  program,  though  it  was  later  simplified  to  remove  the 
quadratic  phase  term  that  appears  in  the  phase  equations.  The  signal  response  phase  equations 
here  shall  be  referred  to  as  the  exact  paint  target  response. 

The  coordinate  system  used  is  shown  in  Fig.  2.1.  The  data  gathering  platform  moves  with 
velocity  V  in  the  -x  direction  and  at  an  altitude  h.  The  ideal  poin»  .arget,  P.  is  located  at 
(x,.y,)  which  is  measured  relative  to  the  reference  point  Q.  located  at  (x„.Y0)2.  Q  is  located  at 
(0.0)  on  the  stationary  axes  x'— y*.  The  terrain  patch  is  assumed  to  be  a  square  of  size  L  by  L. 
Note  that  the  position  of  Q  changes  from  pulse  to  pulse,  since  the  reference  axes  (x-y)  are  mov¬ 
ing  relative  to  Q.  The  range  variables  rtln  and  rlt  are  the  distances  to  Q  and  P.  respectively,  and 
0„  is  the  squint  angle  in  the  ground  plane  during  the  nth  pulse.  The  subscript  n  is  used  to 
describe  variables  which  change  from  pulse  to  pulse  and  so  indicates  the  variable  during  the  nth 
pulse  when  the  look  angle  is  0„.  Later,  n  will  be  dropped,  so  that  each  variable  is  defined  for  an 
arbitrary  angle  0. 


2lt  is  assumed  that  the  distance  from  the  data  plallorm  to  the  reference  Q  is  known  during  the  entire  data  collec 
tion  window.  Errors  in  this  distance  measure  can  delocus  the  reconstructed  image,  but  can  be  corrected  through  a  pro 
cedure  known  as  autolocus  [ihj. 
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Figure  2.1  Coordinate  system  for  Exact  Point  Target  Response  Derivation 


The  transmitted  puise  is  a  high  bandwidth  chirp  signal,  i.e..  a  linearly  modulated  sinusoid: 


s(t) 


COSWtt  +  ^y-)  I  t  I  ^  y 


(2.1) 


0 


otherwise 


where  tue  is  the  carrier  frequency,  v  is  the  linear  FM  sweep  rate  and  T  is  the  puls*  dura* ion. 


$ 


The  nth  received  pulse  is 


8.(0  = 


A  COsfw/t— T„)  +  y(t— Tn)2) 
0 


1 t— r„  1  <T/2 
otherwise 


(2.2) 


where  A  is  the  free  space  propagation  attenuation  factor  and  T„  is  the  two-way  travel  time  of 
the  nth  pulse.  If  Y0»L,  the  propagation  attenuation  of  the  signal  will  be  approximately  the 
same  over  all  the  pulses.  A  will  be  set  to  unity.  Note  the  quadratic  phase  term  in  (2.2).  This 
will  create  problems  later  on.  Note  also  that  this  particular  analysis  assumes  that  the  aircraft 
position  is  essentially  constant  during  each  pulse.  The  analysis  of  Munson  et  al..  [7]  showed 
that  the  aircraft  can  have  an  instantaneous  velocity  of  zero  during  each  pulse,  and  that  (classi¬ 
cal)  Doppler  processing  is  therefore  not  strictly  required  for  SAR. 


Now.  we  apply  stretch  processing,  i.e..  pulse  compression3,  to  the  returned  signal.  This  is 
the  technique  of  compressing  the  chirp  pulse  (and  thereby  increasing  the  range  resolution  of  the 
radar)  by  coherently  mixing  the  return  signal  with  a  simulated  return  from  the  reference  tar¬ 
get.  Typically,  this  signal  is  generated  by  the  transmitter  sweep  oscillator  and  delayed  accord¬ 
ing  to  the  distance  to  the  known  reference.  If  the  distance  to  the  reference  is  not  known  pre¬ 
cisely.  the  simulated  reference  chirp  signal  will  also  be  in  error  and  will  cause  defocusing  of  the 
image,  i.e..  point  targets  will  become  smeared.  The  nth  pulse  reflected  from  the  reference  Q 
will  have  the  form 


gijt)  * 


A  cos<a»t(i~T,Nl)  +  yU-To,,)2) 
0 


1 1— Tm,  I  <Tt/2 
0 


(2.3) 


where  Te  is  the  effective  pulse  width  determined  by  the  size  of  the  terrain  being  imaged. 


Multiplying  g„(t)  and  gou(t).  and  low  pass  filtering  (to  retain  only  the  difference  fre¬ 
quency).  we  obtain  the  phase  function  <t>: 


5  Stretch  processing  and  pulse  compression  have  seemingly  opposite  connotations,  but  refer  to  the  same  technique 
of  processing  a  stretched  pulse  and  compressing  the  pulse-width  through  autocorrelation  of  a  well-designed  pulse  signal, 
e.g..  a  linearly  modulated  chirp  pulse.  See  Klauder  et  al.  (17]  for  more  information  on  chirp  radars. 
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Lowpass  Filter  {  gn(t)  •  gon(t)  }  =  cos(4>n(t»  (2.4) 

where 

4>n(t)  =  <tic(i—Tn)  +  y(t-r0n)2  -  o>c(t-r0n)  -  y(t-r„)2  (2.5a) 

=  (t„~ r0„X(«>e  +  n) -  y (r 2~r02)  (2.5b) 

This  phase  function  corresponds  K  the  phase  of  the  1  'in-phase)  channel  of  the  demodulator 
(mixer  and  low-pass  filter).  To  obtain  the  Q  channel  (quadrature)  output,  multiply  the  return 
signal  by  a  form  of  (2.3)  with  cos  replaced  by  sin  and  then  low  pass  filter  to  retain  the 
difference  frequency. 

1  channel  output:  Real  {e*’’’<t)}  *  cos(4>u(t)) 


Q  channel  output:  Jmag  le^"*0}  =  sin(4>„(t)) 

We  would  like  to  now  convert  the  ensemble  of  return  signals  into  a  two-dimensional 
function  and  show  that  it  is  an  approximation  to  the  phase  of  the  2-D  Fourier  transform  of  a 
point  target.  The  two-way  time  delay  t  to  a  target  at  range  r  is  2r/c.  where  c  is  the  propaga¬ 
tion  velocity  of  the  radar  signal  in  free  space.  Substituting  for  r  in  (2.5b).  the  phase  function 
becomes 


♦„(t)  *  —  (rn-rou)(wc  +  vt)-^(rl?-rt?n) 
c  c2 


(2.6) 


It  is  now  useful  to  make  a  change  in  variables 


*  (tJthJ2  -  1 


(2.7) 


Note  that  for  our  geometries.  Iz„l  «  1.  Applying  the  Pythagorean  theorem  to  Fig.  2.1  and 
making  the  substitution  x„  =  Y«  tanfi,,  we  obtain 

12  (x0  +  x,)2  +  (Y„  +  y,)2  +  h2 


r«« 


x  2  +  Y2  +  h2 


(2.8) 
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r»  F  _  <  ,  2  *n  +  2  Y0  y,  +  x2  +  yt2 

—  I  —It  ■  1  ■  -  ■" ■  ■■■  . . .  ■—  ■ 


_  2  Y„  x,  cot0„  +  2  Y„  yt  +  x»2  +  yt2 


(2.10) 


The  change  of  variables  equations  in  (2.7)  can  be  rearranged  to  obtain 

l“n~run  *  r0n(-v/l  ""+  *n  “  D 


(2.11a) 


rn  ~  r02n  ■  r02  Zn 


(2.11b) 


Since  I  zn  I  «  1.  the  term  -v/T  +  z^  can  be  expanded  into  a  power  series 


/.  ,  .  ,  z,i  z2  z,j 

7i+2n  =  i  +  T  +  T  +  _+  .... 


in  which  only  the  first  2  terms  are  retained.  Thus 


^On 

r»  -  Ton  **  — z — 


(2.12) 


Substituting  (2.12)  into  (2.6)  results  in  an  approximate  phase  function 


♦n<0 


zn  2vr0 

-  |«c  +  vt  - 


(2.13) 


Note  that  r„n  and  z„  are  each  functions  of  n.  as  well  as  0„.  the  angle  at  which  the  nth  pulse  is 
transmitted.  The  notation  can  be  modified  to  express  these  variables  simply  as  functions  of  0. 
This  results  in  the  specification  of  4>n(t)  as  a  function  of  t  and  0. 


A...  r«(0)  z(0)  I  2vr«(0) 

♦(t.0)  *  - ^ -  K»>{  +  VI  -  - - - 


(2.14) 


Now.  4>  is  a  two-dimensional  function  which  can  be  mapped  as  a  function  of  0  and  t,  Note  that 
r0  and  z  are  explicitly  shown  as  functions  of  0.  This  notation  is  simply  to  emphasize  this  func¬ 
tional  relationship. 
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The  variable  t  can  be  functionally  translated  into  a  range  variable  by  defining  R  as  a  func¬ 
tion  of  time 


R(t)  =  Kv  t  +  Kotfsel 


(2.15) 


where  Kv  represents  the  scaling  in  the  radial  (time)  direction  and  Koffset  is  the  offset  from  the 
origin  from  which  valid  icceived  data  begin4..  Solving  (2.15)  for  t  and  substituting  the  result¬ 
ing  expression  into  (2.14)  yields  an  expression  for  <t>  in  terms  of  R  and  9,  which  is  a  two- 
dimensional  function  of  polar  coordinates. 


4>(R.0)  = 


ro(0)z(0) 


*  +  ^ 


R—koffset  2w(l(0) 


A  proper  choice  of  the  scaling  factors  Kv  and  Koffsel  will  simplify  (2.16).  Let 


K 


offset 


K„ 


o><  _  2r0 
v  c 


and  substitute  (2.17)  into  (2.16): 


*(R.e)*^.ro(fl)2(0) 


Substituting  (2.10)  into  (2.18)  gives 


2vY„R 


^(R.O)  *  ~ --(yt  +  x,coi0  +  J^-) 


where  p2  «  x,2  +  yt2.  the  squared  radial  distance  of  P  to  Q. 


(2.16) 


(2.17) 


(2.18) 


(2.19) 


To  eliminate  Kv  and  make  (2.19)  into  a  more  useful  form,  choose  Kv  such  that 


Kv  * 


Yoc 

2  r»  sin© 


Yielding 


(2.20) 


*  The  received  signal  in  (2.2)  is  only  valid  during  a  small  time  interval  which  is  dependent  on  the  nearest  and 
farthest  points  of  the  imaged  terrain.  The  nearest  (and  farthest)  point  is  that  point  which  is  the  minimum  land  max¬ 
imum)  distance  trom  the  data  platform  to  the  terrain  patch  over  the  entire  data  gathering  interval. 
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<t>(R,0)  =  -^(y,sin0  4-  x,cos0  +  ~P^~) 

And  in  rectangular  coordinates  where  U“R  cos{0)  and  V=R  sin(0) 

<KU.V)==-L(xtU  +  ytV  +  £^) 
c*  2yu 


(2.21) 


(2.22) 


This  has  the  same  form  as  the  phase  function  of  the  Fourier  transform  of  an  offset  delta  func¬ 
tion.  The  difference  between  this  expression  and  the  ideal  point  target  response  is  the  addition 
of  an  extra  term  -  the  quadratic  phase  term.  This  quadratic  term  is  dependent  on  p2  which  is  the 
squared  distance  of  the  poi.it  target  to  the  image  center  (reference).  As  the  target  is  moved  out 
radially  from  the  center,  the  quadratic  term  grows  and  causes  a  noticeable  smearing  of  the 
response,  specifically,  a  point  target  widens.  Compare  this  to  the  inverse  Fourier  transform.  F. 
of  a  point  target  at  (x,.y(): 


F(U.V)  *  F-'{8(xt.yt)}  «  J 8(xt.yt)  e^*1' +  yV)  dxdy  «  e**'1’  *  y,v) 


FOJ.V) 


j*  x,li  +  y,v+ 


Jyfi 


**  F(U.V) 


(2.23) 

(2.24a) 

(2.24b) 


where  'V  «  Xr  has  dimensions  radians/ meters2. 
cz 

Recall  that  this  is  only  an  approximation  to  the  exact  point  target  response  since  the  high- 
order  terms  of  z„  were  ignored.  The  exact  point  target  response  can  be  formulated  by  substitut¬ 
ing  (2.11a).  (2.11b).  (2.14).  (2.17).  and  (2.20)  into  (2.6): 


<KR.0) 


2vr,?(6) 


\/l+zU>)- 


(2.25a) 


r<?(0)  «  Y<?  csc20  +  h2 


(2.25b) 
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z(6)  = 


2Y0(x,cotO  +  yt)  +  x,2  +  y,2 


(2.25c) 


And  in  rectangular  coordinates: 


c2  * « 


^l+z(U.V)-r  -z(U.V) 


(2.26a) 


r,?(U.V)  =  Y,r  1  +  4-  h2 


(2.26b) 


z(U.V) 


2Y„(x,^.  +  y,)  +  x,2  +  y,2 


r,2(U.V) 


(2.26c) 


These  equations  relate  the  position  of  a  point  target  at  (x,.y,)  to  its  continuous  complex  phase 
function  4KU.V),  The  data  which  are  obtained  from  the  imaging  system  consist  of  samples  of 
♦  on  a  polar  grid. 

The  phase  function  relates  the  complex  2-B  frequency  pattern  to  the  position  of  the  point 
targets.  This  means  that  a  straightforward  spectral  analysis  of  the  ensemble  phase  function 
(though  sampled)  should  lead  directly  to  the  positions  of  the  point  target.  Superposition  is  also 
applicable  here,  so  that  a  continuous  complex  reflectivity  map  can  be  generated  from  the 
Fourier  transform  or  the  reflected  and  processed  pulses. 


2.2  The  Tomographic  Formulation  of  SAR 


Munson.  O'Brien  and  Jenkins  (7]  have  developed  the  theory  of  SAR  from  a  tomographic 
point  of  view.  They  examine  SAR  as  a  narrow-band  filtered  tomography  problem  using  the 
projection  slice  theorem.  Bernfeld  [18]  later  develops  the  same  sort  of  link  between  the  two 
areas,  although  many  of  the  processing  considerations  seem  to  be  lacking. 

In  the  following  development,  spatial  domain  functions  will  be  lower  case.  Fourier 
domain  functions  will  be  upper  case  and  functions  of  polar  coordinates  will  have  the  subscript 
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p;  i.e..  f(x,y)  is  a  spatial  domain  function  in  rectangular  coordinates  and  Fp(p,0)  is  a  polar 
representation  of  the  Fourier  transform  of  f(x,y).  F(o>x.ci)y)  is  given  by 

on  (jo 

F(ft)x,<«)y)  =  /  / f(x.y )  e"J<x“* +  y“y)dx  dy  (2.27) 

— «>— oo 

and 

oo  oo 

f(x,y)  =  -L  /  /  F(o>x.ojy)  e+liM* +  yu>y)dtox  dwy  (2.28) 

47T  —oo—oo 

A  projection  of  a  multidimensional  signal  is  a  new  function  with  dimension  one  less  than 
the  original  signal.  In  effect,  one  of  the  dimensions  is  integrated  out  via  a  line  integral.  The  pro¬ 
jection  of  f(x.y)  at  an  angle  0  is  given  by 

(JO 

f(u;0)  =  J f(u  cos0  —  v  sin0,u  sin0  +  v  cos0)dv  (2.29) 

— oo 

A  projection  is  a  function  of  one  variable,  u.  with  the  second  parameter,  0.  indicating  the  angle 
of  the  projection.  The  transform  of  p(u:0)  which  is  a  one-dimensional  transform  with  respect 
to  u  is 

OO 

P(w;0)«  /p(u;0)e-lu"du  (2.30) 

— oo 

The  projection  slice  theorem  can  thus  be  succinctly  stated  as 

P(w:0)  **  F,,(w.0)  (2.31) 

The  proof  is  found  in  a  variety  of  sources.  This  particular  statement  of  the  theorem  is  given 
in  [13]. 

» 

The  projection  slice  theorem  says  that  the  Fourier  transform  of  a  projection  of  f(x.y)  at 
an  angle  0  is  equal  to  a  slice  of  F(wx.<ay)  taken  at  the  same  angle  0  and  passing  through  the  ori¬ 
gin.  The  consequence  is  that  the  original  function  can  be  reconstructed  by  filling  in  the  Fourier 
data  space  with  the  transforms  of  the  projection  data  and  then  inverse  transforming  the  result. 
The  entire  Fourier  plane  must  be  constructed  with  an  infinite  number  of  projections,  since  the 
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projection  operation  is  performed  at  discrete  angles. 

In  computed  tomography,  f(x.y)  is  an  attenuation  coefficient,  or  density  function,  which  is 
to  be  constructed  through  non-invasive  means.  Collimated  X-rays  oriented  at  the  projection 
angle  6  perform  the  line  integrals,  i.e..  projections,  of  f(x.y).  Unfortunately,  only  a  finite 
number  of  projections  are  taken,  and  each  projection  is  sampled  so  that  the  resulting  data  set 
lies  on  a  polar  grid  equally  spaced  in  0  and  R.  This  results  in  an  ill-posed  reconstruction  prob¬ 
lem  [19].  Though  the  previous  paragraph  describes  a  method  of  reconstructing  the  original 
function  f(x.y)  from  the  sampled  Fourier  data  (called  Direct  Fourier  Reconstruction),  there  are 
other  techniques  such  as  Convolutional  Back-Projection  (CBP)  and  Algebraic  Reconstruction 
Techniques  (ART)  which  produce  very  good  images.  Current  state-of-the-art  systems  have 
sub-millimeter  resolution.  Convolutional  Back-Projection  has  the  greatest  favor  in  CAT  recon¬ 
struction  because  it  has  the  lowest  computational  burden  for  the  required  image  quality. 

Direct  Fourier  reconstruction  suffers  from  the  same  geometric  difficulties  as  the  spotlight 
SAR  problem.  The  output  data  set  is  presented  on  a  polar  grid  and  must  be  interpolated  to  a 
rectangular  grid  prior  to  the  transformation  step.  Currently  the  computational  requirements 
and  resulting  images  are  not  competitive  with  CBP. 

In  the  Munson-O'Brien-Jenkins  formulation  of  SAR[7],  the  approach  to  the  analysis  is 
different  from  that  of  Weis[l5].  It  is  simplified  because  the  altitude  of  the  aircraft  is  zero 
(though  this  is  corrected  later  in  the  paper),  and  the  ground  patch  is  circular  (leading  to  a 
simplified  expression  for  determining  what  part  of  the  return  signal  is  valid).  The 
simplifications,  however,  lead  to  great  insight  into  the  imaging  process.  The  following  analysis 
is  taken  from  [7]  and  the  geometry  of  the  problem  is  in  reference  to  Fig.  2.2.  in  which  the 
radar  is  imaging  a  circular  patch  of  radius  I.,  at  a  distance  R  from  the  patch  center.  The  signal 
transmitted  is  as  before  in  (2.1): 


0 


otherwise 
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2R  2R 

r<>  =  A  I  g(x„.Y0)  I  -cos  ^(t-— -)  +  y(t— — r^-)2  +  arg{g(x0.Y0)} 

V  4*  V 


(2.34) 


Again,  just  as  in  the  previous  analysis,  the  return  signal  contains  a  quadratic  phase  term.,  The 
attenuation  constant  A  can  again  be  set  to  unity  since  it  will  be  approximately  the  same  for  all 
pulses  if  R»l...  Equation  (2.34)  can  be  simplified  to 


r„  *  Re 


2R 

jg(  x». Y0)'s(  t—-~  )dx  dy 


(2.35) 


This  is  the  return  from  a  small  differential  element  at  a  distance  R„  from  the  transmitter.  For 
the  geometry  of  Fig.  2.2,  the  locus  of  points  at  a  distance  Rl(  from  the  data  platform  is  an  arc 
centered  at  the  radar  system.  The  return  signal  corresponding  to  the  superposition  of  all 
differential  patches  at  distance  R„  is  the  line  integral  along  the  same  arc.  However,  if  R»L, 
then  the  circular  (spherical)  radar  waves  can  be  approximated  by  a  plane  wave  normal  to  the  u 
axis  (the  axis  of  transmission  for  an  angle  0).  This  line  integral  is  the  value  of  the  projection  of 
g(x.y)  onto  the  u-axis  at  u„.  Equation  2.35  can  be  thus  rewritten  as 


r„(t) »  Re 


PgU,)’S 


t  - 


2(R  +  U|i) 


Wu 


(2.36) 


where  R«,  has  been  rewritten  as  R.  the  distance  to  the  terrain  center  plus  u«.  Note  that  these 
quantities  are  dependent  on  9  since  R  is  changing  as  the  plane  flies  past  the  patch.  The  return 
signal  from  the  entire  terrain  patch  is  given  by  integrating  r0(t)  over  the  entire  area  (here  the 
assumption  that  the  terrain  Is  circular  simplifies  the  integration  a  great  deal,  though  we  could 
assume  that  a  square  patch  is  circumscribed  by  the  circle,  and  the  terrain  between  the  circle  and 
square  has  reflectivity  zero.) 


r‘#(t)  =  Re 


fptAu)'. 


It  - 


2(R  +  u) 


du 


(2.37) 
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Because  the  transmitted  signal  is  a  chirp  pulse,  stretch  processing  is  applied  as  before, 
which  essentially  translates  positional  information  into  frequency  information.  This  is  done  by 
mixing  (multiplying)  this  return  r'#  signal  with  a  reference  chirp  from  the  ground  patch  center 


cos[wr(t— r„)  +  y(t— tJ2] 


(2.38) 


2R 

where  r«  *  — the  two-way  travel  time  to  the  terrain  center.  This  is  low-pass  filtered  to 
obtain 


r<f(t)  *  yRe 


f  pa(u)‘¥(i)-Y(t)du 


(2.39) 


where 


♦<i)«e 


(2.40a) 


is  the  quadratic  phase  term,  and 


Y(t)  *  exp 


— j— (o»t  +  ^-(t— r0))u 

C  2 


(2.40b) 


The  quadrature  component  is  obtained  by  multiplying  by  the  phase  shifted  reference  chirp 


sin[<Dc(t-T„)  +  y(t-To)2] 


(2.41) 


and  low  pass  filter,  giving  the  imaginary  portion 


r^(t)  *  ylm  I f  pg(u)’¥(i)’Y(i)du 


(2.42) 


of  a  complex  signal 


Cod)  =  i 


/pe(u)-¥(t)-Y(t)du 


(2.43) 
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If  the  quadratic  phase  term,  ¥(t).  can  be  removed,  C#(t)  resembles  the  Fourier  transform  of 
Po(l)  in  the  u  coordinate  plane,  evaluated  at  [toc  +  y(t— r0)]. 

Cj(t)  =  1p<K  +  y(t-To)]  (2.44) 

If  the  transmitted  chirp  pulse  had  an  infinite  bandwidth,  then  C'#(t)  would  represent  a 
continuous  sample  (slice)  of  the  Fourier  transform  of  the  ground  patch  reflectivity,  g(x.y). 
(Note  that  if  the  terrain  reflectivity  is  isotropic  i.e..  g(x.y)  is  not  a  function  of  the  radar  look 
angle,  then  the  projection  is  symmetric  about  the  origin.)  However,  the  pulse  is  only  of  finite 
width,  thus  finite  bandwidth,  and  so  the  processed  return  represents  only  a  small  section  of  the 
Fourier  slice.  The  return  pulse  is  only  valid  for  the  interval 

asm-T  **-';> +  T  a«) 

C  2  c  2 

This  interval  represents  the  two-way  time  for  the  beginning  of  the  pulse  to  travel  to  the  god  of 
the  terrain  patch  until  the  gnd  of  the  transmitted  pulse  reaches  the  beginning  of  the  terrain 
patch.  These  calculations  are  based  on  minimum  and  maximum  distances,  which,  for  a  rec¬ 
tangular  patch,  change  from  pulse  to  pulse  (L  would  change  non-linearly  with  0).  Again,  the 
round  patch  has  simplified  the  calculations.  Substituting  l  of  (2.45)  into  (2.44)  yields  an  ine¬ 
quality  for  the  range  of  the  argument  of  P* 

i(«t-yT+^)  <  K  <  i(W<  +  vT-l^)  (2.46) 

c  c  c  c 

For  a  typical  transmitted  chirp  pulse  where  w{±sT»2sl./c.  (2.44)  reduces  to 

l(wt-ij.)  <  Ru  <  i(wc  +  il)  (2.47) 

c  2  c  2 

Rw  is  bounded  by  the  inner  and  outer  radii  of  the  slice  of  the  Fourier  transform  of  the  terrain 
reflectivity  represented  by  the  processed  return  signal.. 

Thus,  as  before,  the  coherently  processed  signal  is  a  section  of  a  slice  of  the  Fourier 
transform.  This  signal  is  also  sampled  uniformly  in  range  so  that  the  collection  of  sampled 
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returns  assumes  a  polar  grid,  offset  from  the  origin  and  having  an  angular  sweep  width  of  0M. 

Many  issues  not  covered,  such  as  the  quadratic  phase  term,  modification  for  non-zero  plat¬ 
form  height,  curvature  of  the  radar  signal  and  the  effects  of  plane  movement  (inducing  a 
Doppler  term)  are  not  included  here,  but  can  be  found  in[7].  As  before,  the  quadratic  phase 
term  has  the  effect  of  smearing  targets  which  are  distant  from  the  patch  center..  This  is  obvi¬ 
ous.  since  the  u2  term  represents  the  square  of  the  radial  distance  of  a  target  from  the  terrain 
origin.  It  is  the  same  quadratic  phase  term  found  in  the  Weis  analysis  [15]  and  is  difficult  to 
remove  during  processing.  It  can  be  shown  that  this  term  limits  system  resolution  and  terrain 
patch  size. 

23  Comparison  of  CAT,  SAR,  and  Analyses 

Munson  el  al.  [7]  compare  CAT  and  SAR  and  note  important  differences: 

( 1 )  As  a  result  of  stretch  processing,  the  processed  SAR  signal  is  the  Fourier  transform  of  the 
projection.  In  CAT.  the  projection  is  still  part  of  the  spatial  domain. 

(2)  The  processed  SAR  signal  gives  only  a  small  part  of  the  projection  transform  devoid  of  dc 
components,  which  can  result  in  the  loss  of  edges  and  sections  of  constant  reflectivity. 

(3)  In  SAR.  the  projection  is  taken  normal  to  the  axis  of  imaging  rather  than  parallel  to  the 
axis  (as  in  CAT). 

These  two  different  ways  of  looking  at  the  SAR  problem  lead  to  essentially  the  same 
result:  the  collected  data  are  approximately  samples  of  the  Fourier  transform  of  the  terrain 
reflectivity  function.  With  Doppler  interpretation  (Weis),  the  image  is  obtained  by  spectral 
estimation,  i.e..  a  forward  FFT.  to  translate  sinusoids  into  discrete  range  bins.  In  the  tomo¬ 
graphic  formulation,  the  data  represent  samples  of  the  Fourier  transform  of  the  complex 
reflectivity  which  must  be  inverse  transformed  to  obtain  samples  of  the  complex  reflectivity.. 
Both  interpretations  are  correct,  but  can  lead  to  confusion  about  the  role  of  Doppler  processing 
and  the  direction  of  the  Fourier  transform.  In  practice,  it  is  of  little  concern,  since  the  phase  is 
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discarded  in  the  final  displayed  image  and  forward  or  inverse  transforms  result  in  the  same 
(though  inverted  and  scale)  magnitude  map. 

It  should  be  pointed  out  that  the  Direct  Fourier  Reconstruction  approach  is  but  one  algo¬ 
rithm  for  the  terrain  mapping.  Farlier  versions  of  strip  mapped  SAR.  based  on  the  same  princi¬ 
ple  used  the  concept  of  multiple  correlators  to  perform  the  stretch  processing.  The  correlators 
performed  the  job  of  matched  filters,  which  were  tuned  to  point  targets  in  each  range  bin.  hav¬ 
ing  the  impulse  response,  which  was  the  complex  conjugate  of  the  response  of  a  point  target  of 
that  range  bin.  This,  however,  invojves  much  more  processing  than  the  FFT  approach,  though 
hardware  implementations  have  been  designed  and  built  which  can  reduce  the  processing  lime 
through  a  pipelined  architecture  (20]. 

2v4  Other  Inversion  Techniques 

The  two-dimensional  interpolation  scheme  is  the  most  direct  method  for  the  SAR  image 
reconstruction  problem.  Since  a  relationship  was  formed  between  SAR  and  CAT.  it  can  be 
theorized  that  some  of  the  tomographic  reconstruction  algorithms  may  be  used  in  SAR.  Unfor¬ 
tunately.  many  of  them  depend  on  the  full  set  of  projection  data  (0M  ■  2ir)  and  thus  are 
difficult  to  apply  to  SAR.  In  particular,  the  algorithms  that  may  be  of  use  are  convolutional 
back-projection  (CBP)  and  the  Hankel  transform. 

2A1  Convolution  back  projection  and  the  Hankel  transform 

Quite  a  bit  of  research  has  been  done  in  the  area  of  tomographic  reconstruction  algorithms 
and  their  computational  burden  versus  image  quality  (similar  to  the  work  presented  here  for 
SAR).  A  celebrated  paper  by  Pan  and  Kak  (21}  details  the  tradeoffs  of  direct  Fourier  inversion 
with  interpolation  and  CBP.  but  only  very  simple  intei  polators  are  used,  notably  the  nearest 
neighbor  and  bilinear  (used  with  a  modification  of  the  zero-padding-FFT  technique  described  in 
Chapter  3).  Interpolation  via  the  circular  sampling  theorem  (see  Stark  (5. 22]  and  Fan  [23]  for 
further  discussion  about  this  theorem)  is  presented,  but  again,  this  requires  a  periodic  data  set 
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(complete  circular  arcs)  to  be  useful.  Heffernan  and  Bates  [24]  also  compared  the  very  low- 
order  interpolator  methods,  but  applied  them  to  the  projection  data  prior  to  back-projection, 
rather  than  2-D  direct  inversion.  Similar  work  has  been  done  by  Mersereau  [13].  Desai  and 
Jenkins  [25.  26]  have  researched  the  use  of  the  CBP  algorithm  for  the  SAR  problem  and  have 
achieved  image  quality  comparable  to  the  direct  inversion  technique,  although  at  the  expense  of 
additional  computation  order.  Even  if  the  SAR  problem  could  be  considered  a  limited  view 
tomography  problem,  the  reconstruction  techniques  of  bandlimited  extrapolation  to  supply  the 
rest  of  the  Fourier  data  set  would  be  far  too  costly  for  real  time  reconstruction. 

The  Ilankel  transform  allows  for  direct  polar-to-polar  image  reconstruction  from  the 
sampled  Fourier  data  by  expanding  the  two-dimensional  image  into  a  Fourier  series  and  then 
performing  a  Hankel  transform  on  these  coefficients  to  obtain  the  Fourier  series  coefficients  of 
the  image.  The  use  of  the  Hankel  transform  to  invert  the  polar  Fourier  data  set  directly  has 
been  met  with  limited  success  [27]  due  to  the  cost  of  computing  the  needed  transforms.  Work 
is  being  done  to  develop  a  Fast  Hankel  Transform  [2S].  but  this  transform  still  requires  a  full 
circular  data  set. 

2AJ2  Spectral  analysis  of  nonuniformly  sampled  data. 

Hostetter  [29.  30]  developed  a  control  systems  approach  to  spectral  observation  of  data 
which  is  irregularly  sampled.  The  algorithm  is  not  practical  for  large  data  records,  however, 
because  it  is  of  order  N2  for  N  input  points.  For  an  N  by  N.  2-D  array,  this  would  require 
(XN4)  calculations,  which  is  not  competitive  with  even  the  most  sophisticated  2-D  interpolator. 

2.5  Polar-Rectangular  Geometry  for  Interpolation 

The  first  two  sections  of  this  chapter  demonstrated  how  the  spotlight  mode  SAR  data  set 
falls  on  a  polar  grid  of  limited  size.  The  major  thrust  of  the  work  here  is  to  use  the  Direct 
Fourier  Reconstruction  algorithm  to  generate  the  complex  reflectivity  map  of  the  imaged  terrain 
and  to  study  effects  of  different  two-dimensional  interpolators  available.  Fig.  2.3  details  the 
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polar-to-rectangular  grid  relationships  and  the  parameters  needed  to  specify  the  relative 
geometry  between  the  two  rasters.  The  input  polar  grid  is  specified  by 

(1)  Rmm  the  inside  radius  of  the  torus  of  known  Fourier  data. 

(2)  the  maximum  look  angle  of  the  SAR  system. 


Figure  2.3  Polar-lo-Rectangular  Grid  Geometries. 
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(3)  Ar,  the  radial  sampling  increment. 

(4)  A0  the  angular  sampling  increment. 

The  output  rectangular  (square)  raster  is  specified  by 

( 1 )  uO.vO,  the  bottom-left  coordinate  of  the  square  grid. 

(2)  Au.Av  the  sampling  increments  for  each  of  the  u  and  v  dimensions.  Typically,  these  will 
be  equal. 

The  resolution  of  the  system  is  governed  by  the  width  of  the  Fourier  section  in  each 
dimension.  These  widths,  in  turn,  are  determined  by  the  radar  carrier  frequency  wc.  chirp 
sweep  width,  and  the  maximum  look  angle.  0M.  In  the  range  direction,  the  width  is 
K,nax  ”  R|im>*  the  bounds  on  the  range  variable  given  above  determined  by  the  radar  chirp  sweep 

2vT 

width.  From  (2.47)  above,  the  resolution  is  given  as  approximately  — - —  In  the  azimuthal 

dimension,  the  v-width  is  approximately  that  of  the  azimuthal  extent  of  the  annuluar  section, 
approximately  4o»{  sin0M.  If  the  processed  response  signal  represented  the  Fourier  space  exactly 
(no  quadratic  phase  term),  then  a  point  target  placed  at  (x„.Y0)  would  yield  a  two-dimensional 
sine  due  to  the  limited  extent  of  the  Fourier  space  available.  The  resolution  of  the  system  can 
be  defined  as  the  width  of  the  sine  mainlobe  divided  by  2  since,  theoretically,  this  represents 
the  separation  of  two  resolvable  point  targets.  In  practice,  the  resolution  is  degraded  though 
constructive  and  destructive  interference  of  the  phase  function. 

To  achieve  maximum  resolution,  it  is  important  to  interpolate  from  the  polar  grid  to  the 
largest  square  which  inscribes  the  angular  section  (geometry  shown  in  Fig.  2.4).  The  range 
resolution  is  determined  by  the  chirp  sweep  width  (time-bandwidth  product),  and  the  azimu¬ 
thal  resolution  is  determined  by  the  center  frequency  otc  and  the  look  angle  0M-  These  parame¬ 
ters  specify  the  annular  section  location  and  size  in  the  Fourier  domain.  The  parameters  used  in 
the  simulations  were  wc  and  the  look  angle  0M,  from  which  the  chirp  sweep  width  and  square's 
size  and  position  are  determined.  The  solution  in  this  case  is  straightforward  trigonometry. 


I^et  W  be  the  width  of  the  square  to  be  determined.  Rmill  and  Rmax  are  also  to  be  found  -  relat¬ 
ing  back  the  chirp  sweep  width  of  the  radar. 


Figure  2.4  Determining  the  Largest  Square  Constrained  by  <uc  and 


«*  vireni 
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There  is  a  very  slight  gap  between  the  right  edge  of  the  square  and  the  torus  on  the  v«0  axis. 
Rmax  is  given  by 


Rmax  *“ 


|(rmi„  +  W)2  +  (~)2 


since  the  square  must  touch  the  outer  radii  for  maximal  size.  W  is  related  to  Rmin  by 

9M 

W  =  2Rmilltan~ 

Substituting  (2.50)  into  (2.49)  gives 

R„2  = 

But  from  (2.43).  RmM  «  2q0  —  Rmill.  Substituting  this  into  (2.51 )  results  in 


d  6 

'max  ®  ( Rmln  4"  2Rm,nlan  — — —  )2  +  (Rmmtan-y  )2. 


0  0 

(2q<i  "*  Rmin^'  =  (Rmin  2Rmmtan-y  )“  +  (Rmii>tl*n-~~)2  . 

Finally,  solving  for  Rmju  yields  a  quadratic  equation  which  has  two  roots: 


4q0±4q„ 


R 


mm 


(1  +  2tan-^  )2  +  tan2~ 
2  2 


2(1—0  +  2tan-Jl.)2  +  tan2-il) 
2  2 


(2.49) 


(2.50) 


(2.51) 


(2.52) 


(2.53) 


The  positive  solution  yields  the  correct  result  and  Rm#x  and  W  can  be  obtained  via  (2.48)  and 
(2.50). 

In  this  way.  the  sweep  width  of  the  radar  system  can  be  determined  for  a  given  look  angle 
(assuming  the  center  frequency  remains  constant).  If  the  chirp  bandwidth  is  less  than  that 
given  above,  the  range  resolution  will  deteriorate:  if  it  is  more,  then  the  extra  data  are  not  used. 


2.6  Oversampling  and  Presumming 


A  typical  SAR  system  samples  the  demodulated  output  at  a  much  higher  rate  than  would 
be  necessary,  based  on  the  Nyquist  rate  for  the  recorded  signal.  This  is  because  the  antenna 
pattern  is  not  an  ideal  step  function  which  is  limited  to  the  terrain  being  imaged.  Rather,  it  has 
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.  transition  band  and  sidelobmt  which  can  detect  targets  outside  the  desired  patch.  These  tar- 
gt  *  are  al.ased  back  into  the  processed  image  and  result  in  false  targets.  By  sampling  at  a 
***«  rate,  and  then  digitally  hit,™,  the  data  prior  intention,  these  outside  tar*,*  can 

^  SU“'  SamP’i”S  lh'  — *  done  for  the  same  r«son.  tho„sh  it  is 
referred  to  as  Dopp.e,  oversampling.  si„ct  ,he  azimuth  antenna  pattern  is  obtained  through 

Doppler  methods  (31).  The  oversampled  range  lines  are  then  low-pass  titered  in  the  azimuth 
direction  narrow  ,h,  antenna  beamwidth  and  a, so  to  reduce  ,h,  system  noise  , coherent 
summing  sever.,  range  lines  into  one,.  „  Wou,d  smm  that  th,  pretitered.  highiy  overtamp, ed 
data  would  be  better  input  the  interpolation  stag,.  si„ce  the  Coser  spaced  dau  would  mduce 
•he  interpolation  error.  This  is  no,  don,  b«.»se  the  storage  requirements  for  such  vohnnes  of 
dau  are  prohibitive.  This  topic  is  treated  by  Brown  *  of.  i„  designing  „ptimal  presummj„g 

dUers  which  ret.m  ezimutha,  reso.utKm  and  a„ow  for  limited  dam  storage  (32).  Haymr  (33) 

U*s  the  azimuthal  oversampling  on  a  keystone  grid  to  analyse  the  noise  elfects  of  a  one- 

dunensiona,  nearest  neighbor  interpolator  and  to  pmpos,  an  adaptive  presumm.ng  fdter  in  the 
azimuthal  direction. 


THE  INTERPOLATION  PROBLEM 


Before  discussing  two-dimensional  interpolators,  it  is  useful  to  examine  the  one¬ 
dimensional  case  and  see  what  can  be  extracted  from  there.  Past  experience  has  shown,  how¬ 
ever.  that  generalizing  from  one  to  multiple  dimensions  in  signal  processing  algorithms  is  not  as 
simple  as  adding  a  subscript.  The  problem  worsens  when  the  input  grid  is  non-rectangular.  and 
becomes  still  worse  when  it  is  irregular  (the  sample  spacing  cannot  be  parametized).  First,  we 
will  look  at  what  approximation  theory  can  give  us  and  then  some  classical  interpolation  con¬ 
cepts.  Local  versus  global  interpolation  is  discussed,  as  well  as  separable  versus  non-separable 
interpolation.  Then  the  DSP  approach  to  interpolation  will  be  reviewed.  Most  classical 
approximalion/interpolation  theory  texts  deal  with  real  data  only,  although  most  algorithms 
can  be  expanded  to  the  complex  domain  without  loss  of  functionality. 

3.1  Classical  Interpolation 

Classical  approximatum  theory  attempts  to  solve  the  following  problem  [34] : 

If  V  is  a  normed  linear  space  and  W  is  a  subset  of  V,  then  given  a  v  C  V.  find  a 

w’  C  W  such  that 

II  v  —  w*  II  ^  II  v  —  w  II 

for  all  w  C  W.  where  II  •  I!  is  a  linear  norm  on  the  space  V 

The  norm  is  typically  the  Chebyshev  norm  (f  or  the  Euclidian  norm  (L2).  In  our  case,  V  is 
the  set  of  two  dimensional  functions,  and  W  is  the  set  of  spatially  limited  functions,  e.g.,  the 
ideal  low-pass  filter  (bandlimited)  is  approximated  with  a  truncated  or  tapered  sine  function 
(spatially  limited,  but  not  bandlimited  due  to  the  truncation). 

Classical  interpolation  theory  sets  out  to  solve  the  following  problem: 


Given  an  m  partition  on  tbe  interval  l4a.b]  with  partition  set  X  =  Xj.x2.  *  •  •  .xm 
such  that  a=x1<x2<  •  •  •  <xm=b.  Let  y,  =  f(xj).  Find  a  function  g*(x)  such  that 

g*(xj)  =  yj  l^i^m 

and 

II  f(x)  —  g*(x)  II  <  II  f(x)  —  g<x>  If  for  Xj  <  x  <  xN  . 

This  allows  for  a  very  general  specification  of  g(x)  which  gives  rise  to  the  use  of  piecewise 
functions  in  representing  g(x).  i.e..  piecewise  polynomials  and  splines3.  If  the  subspace  W  and 
V  are  the  same,  then  an  exact  interpolanl  can  be  found.  i.e..  the  function  can  be  reproduced 
exactly.: 

It  is  worth  noting  a  few  points  about  approximation  and  interpolation  for  DSP  applica¬ 
tions.  With  interpolation  theory,  we  are  required  to  obtain  a  function  which  passes  through 
each  of  the  original  datum.  This  has  the  disadvantage  that  a  noise  corrupted  signal  will  be 
reconstructed  with  a  function  passing  through  each  noisy  sample,  rather  than  smoothing  it  out. 
If  the  condition  that  g*(t)=yi  is  relaxed,  then  we  can  produce  an  approximating  function  g*(t) 
which  minimizes  the  mean  squared  error  between  g*(t)  and  the  known  data  points.  Let  the  sam¬ 
ple  error  e(  be  defined  as 

«.  =  y,-  g(x.) 

Find  a  function  g*(t)  which  minimizes 

m 

EtfO-E'e.12 

1=1 

The  minimizing  function  g'(t)  is  our  interpolant.  In  this  case,  however,  g’(x,)  ^  y,  in  general, 
and  the  original  data  will  not  be  reproduced  by  g, 

5  Throughout  this  chapter,,  the  interpolating  functions  will  be  written  with  its  argument  as  t  or  x  interchangeably 
because  x  is  borrowed  1  rom  the  mathematics  area  and  t  from  the  engineering  area. 
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For  approximating  a  continuous  function  (signal),  we  are  usually  required  to  approximate 
it  only  within  a  finite  interval.  This  means  that  our  data  record  is  necessarily  a  fixed  length 
record.  i.e..  it  is  time  (spatially)  limited.  If  our  sequence  is  of  infinite  length,  as  in  sample  rate 
changing  a  voice  or  communications  signal,  then  it  must  be  processed  in  blocks.  Also,  for  a 
fixed  length  record  (such  as  our  SAR  problem)  and  without  the  use  of  signal  extrapolation,  the 
interpolated  reconst  ruction  will  be  only  an  approximation  to  the  sampled  bandlimited  function. 
The  problem  of  processing  a  time-limited  signal  as  if  it  were  also  bandlimited  is  ever  present. 
Recall,  too.  that  these  interpolation  concepts  are  presented  in  the  time  or  spatial  domain,  but  are 
actually  applied  in  the  SAR  frequency  domain  where  duality  is  used  to  justify  assumptions 
about  the  signal. 

3.1.1  One-dimensional  DSP  interpolators 

The  following  section  deals  with  interpolation  from  a  DSP  approach  [35].  It  demonstrates 
that  interpolation  (by  a  rational  factor)  is  really  a  filtering  operation.  The  analysis  is  also  car¬ 
ried  out  assuming  that  the  interpolation  is  done  in  the  time  domain. 

Assume  that  a  continuous  time  function  f(t)  has  finite  energy  and  is  cr  bandlimited.  i.e.. 

£  I  f(t)  1 2  <  oo  (3.1a) 

t*— oo 

and 

F(w)  =  0  I  co  I  >cr  .  (3.1b) 

where  F(g>)  represents  the  usual  Fourier  transform  of  the  function  f(t): 


F(w)  =  J  f(t)-e-J“‘dt 

— oo 

If  f(l)  is  sampled  to  produce  the  infinite  sequence  x(n)  (n,  an  integer), 

x(n)  =  f(nT)  —oo  <  n  <  oo 


(3.2) 


(3.3) 


where  T  is  the  sample  period  and  T  <  tr/cr,  then  the  function  f(t)  can  be  exactly  reproduced 
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with  the  infinite  sum 


f(t)=  Z  x(k). 


sin 

-J-(t-kT) 

T 

k*-~  y(t-kT) 


(3.4) 


Equation  (3.4)  can  be  recognized  as  the  classical  reconstruction  theorem  for  bandlimited  signals. 

It  is  customary  to  refer  to  the  interpolation  kernel  in  (3.4)  as  the  function  sinc(r)  =  S'n-T -  so 

T 

that  (3.4)  becomes 


f(t)  =  Z  x(k)sinc(X(t-kT)).  (3.5) 

k  *  -on  I 

The  reconstruction  described  by  (3.5)  is  equivalent  to  passing  the  impulse  sequence  x(n) 
through  an  ideal  low-pass  filler  with  cutoff  frequency  a>  *  rr/T.  The  low-pass  filler  removes 
the  copies  (aliases)  of  F(«)  which  are  replicated  every  2w7T  by  the  sampling  operation. 

Schafer  and  Rabiner  [36]  showed  that  in  the  realm  of  digital  signal  processing,  interpola¬ 
tion  is  a  type  of  filtering  operation.  Interpolating  by  an  integer  factor  L*T/T  is  accomplished 
through  first  forming  a  new  sequence.  v(n).  by  inserting  L-l  zeros  between  the  input  data 
points.  x(n)  from  (3.3).  and  then  filtering  this  resulting  sequence  with  an  ideal  low-pass  (LP) 
filter  (appropriately  scaled). 


v(n)  = 


x(n/L) 

0 


n  *  0,±L.±2L.  •  •  • 
otherwise 


(3.6) 


y(n)=  Z  v(k)-h(n-k)  (3.7) 

k=-v» 

where  h(n)  is  the  impulse  response  of  an  ideal  low-pass  filter. 

The  output  sequence  y(n)  will  be  exact  samples  of  the  original  function  f(t)  sampled  at  T 
where  T'  =  T/L.  When  the  input  sampling  rate  T  is  greater  than  the  Nyquist  rate  n/<r,  the 
requirements  of  the  LP  filter  are  relaxed  so  that  all  that  is  required  is  an  LP  filter  in  which  the 
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passband  is  approximately  unity  for  rr/ T  <  1  o>  I  and  approximately  0  for  lot  I  >o\  Optimal 
FIR  filters  can  easily  be  designed  which  meet  these  criteria  [37].  however,  there  may  be  an 
additional  requirement  that  the  original  samplp  values  remain  undisturbed  by  the  filtering 
operation.  This  additional  constraint  may  degrade  the  filter  a  bit  by  causing  less  attenuation  in 
the  stopband  [35]. 

If  the  input  sequence  is  sampled  close  to  the  Nyquist  rate,  the  LP  filter  must  then  be  very 
close  to  the  ideal  in  order  to  avoid  aliasing.  This  usually  necessitates  a  high-order  filter  which  is 
very  costly.  This  has  the  same  effect  as  applying  (3.5)  and  then  sampling  the  continuous  out¬ 
put  f(t)  with  a  sampling  period  T*  (although  the  continuous  function  f(t)  is  of  course  never 
really  formed  in  a  finite  state  digital  computer).  When  T/T'  *  L/M  is  rational,  then  we  can 
convert  to  the  new  sample  rate  by  first  interpolating  by  a  factor  of  L.  low-pass  filtering,  then 
decimating  by  M  (retaining  only  every  Mth  point). 

A  difficult  problem  occurs  when  T/T*  is  irrational.  In  this  case,  we  may  never  have  a 
point  at  which  x(n)  and  y(m)  coincide,  and  so  the  above  algorithm  cannot  be  implemented.  We 
must  resort  to  (3.4)  for  taeh  point  y(m)  (see  Fig.  3.1).  The  infinite  sum.  however,  is  impracti¬ 
cal.  and  thus  we  form  a  suboplima!  interpolation  kernel  h(t)  which  has  finite  support. 

y(n)  *  f(nT)  -  £  x(k)  h(nT-kT)  (3.8) 

k»-K 

The  notation  f  indicates  that  the  reconstruction  is  approximate.  In  the  more  general  case  where 
y(0)  occurs  at  some  toffst,  such  that  y(O)  *  f(w„,).  then  (3.8)  becomes 

y(n)  *  fUtf^+nr)  (3.9a) 

K 

=  £  x(k)h(toffsc,+nT-kT)  (3.9b) 

k=-K 


This  suboptimal  filter  may  not  remove  all  of  the  energy  at  frequencies  above  cr,  and  will  conse¬ 
quently  generate  errors  in  the  continuous  reconstruction  f(t),  and  so,  in  the  sequence  y(m). 
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tutor  , 

x(e)  x6) 

*(i) 

x(h)  x(0 

output 

ft*)  y(») 

y(3> 

v(h)  y(c)  y(*>)  y(*) 

. . * - *  " 

Figure  3.1  Irrational  Rate  Change  with  Output  Offset. 

In  the  simpler  integer  ratio  rate  conversion.  h()  is  an  FIR  filter  specified  in  the  digital 
domain  as  a  sequence.  Here,  the  function  h(t)  must  be  a  continuous  function  defined  within  the 
limits  of  the  summation,  rather  than  a  discrete  sequence.  This  is  a  consequence  of  convolving  a 
sampled  sequence  with  a  continuous  function  when  the  output  samples  do  not  match  up  with 
the  input  sample  spacings  and  the  rate  change  is  irrational.  The  continuous  time  filter  must 
have  a  transform.  H(o»).  which  removes  the  periodic  copies  of  F(o»)  which  are  above  ir/T.  The 
performance  of  the  interpolator  Is  determined  by  how  closely  the  Fourier  transform  of  h(t) 
approximates  the  ideal  LP  filler.  That  part  of  the  spectrum  which  is  not  removed  above  ir/T 
will  not  simply  show  up  as  high-frequency  noise,  but  will  be  aliased  back  into  the  baseband  in 
the  spatial  domain  when  this  continuous  reconstruction  is  resampled.  This  will  be  discussed 
more  in  Chapter  4. 
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In  theory,  an  arbitrary  conversion  rate  R  could  be  approximated  by  a  rational  rate  change 
L/M,  but  as  L  and  M  become  large  so  as  to  approximate  an  arbitrary  rate  change,  the  fillers 
grow  to  a  very  high-order.  As  an  example,  a  rate  change  of  1.045  requires  L/M  -  209/200.  that 
is.  interpolating  by  a  factor  of  209  and  then  retaining  only  every  200lh  point.  At  the  expense 
of  sample  delay,  processing  can  be  broken  into  multiple  stages,  each  of  smaller  order.  Several 
methods  have  been  designed  for  breaking  down  a  high-order  LP  filter  into  several  low-order 
stages  where  the  decimation  is  reduced  by  factors  of  two  [38.  39],  or  reduced  to  two  stages, 
each  of  lowered  order{40].  Rabiner  and  Crochiere  present  an  optimization  technique  for  break¬ 
ing  down  the  L/M  inierpolator/decimator  (I-D)  into  a  series  of  K  smaller  order  FIR  filters 
where  the  optimization  is  done  to  minimize  the  total  number  of  multiplies  and  adds  per  second 
(MADS)  [35. 41].  They  have  shown  that  by  breaking  down  a  high-order,  single-stage  L/M 
interpolator/decimator  0-D)  into  several  small  order  L/M  1-D  sections,  a  significant  computa¬ 
tional  savings  can  be  achieved.  This  notion  can  be  generalized  to  design  more  efficient  low-pass 
filters  through  several  interpolation  and  decimation  stages  [42].  It  should  be  cautioned,  how¬ 
ever.  that  cases  can  arise  in  which  there  is  no  savings  in  breaking  down  the  single  stage  [43]. 

Ramsiad  [44]  has  developed  structures  for  conversion  between  arbitrary  sample  rates.  He 
discusses  a  hardware  structure  which  does  real  time  calculation  of  the  interpolator  coefficients, 
which  is  essentially  a  lime  varying  filter.  He  also  presents  a  "hybrid'’  interpolator  that  does 
rational  sampling  rate  conversion,  as  previously  discussed,  followed  by  a  nearest  neighbor, 
linear,  or  possibly  LaGrange  interpolator,  to  generate  data  which  fall  between  those  equi-spaced 
output  samples.  The  approach  Is  hardware  oriented  with  very  specialized  structures  to  gen¬ 
erate  the  lime-varying  coefficients.  While  this  may  be  less  of  i  concern  for  current  VSI.l  tech¬ 
nology.  it  would  be  useful  to  have  an  algorithm  that  could  be  implemented  in  a  readily  avail¬ 
able  serial  or  parallel  processor  or  one  of  the  commercially  available  DSP  chips. 
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3.2  Global  versus  Local  Interpolators 

Nearly  all  interpolators  that  are  used  have  local  support.  That  is.  the  interpolation  ker¬ 
nels  are  finite  in  length  and  use  input  data  in  only  a  small  neighborhood  of  the  output  point. 
The  advantages  to  this  class  of  interpolators  are  high  speed  and  low  cost.  They  are  faster  and 
easier  to  implement  because  only  a  small,  limited  number  of  data  points  are  involved  in  the 
summation  (convolution).  The  disadvantage,  of  course,  is  a  poor  signal  reconstruction. 

From  the  previous  sections,  we  know  that  the  ideal  interpolator  impulse  response  is  the 
sine  function  which  is  of  infinite  length.  This  is  called  a  global  interpolator  because  it  uses  sam¬ 
ple  data  (albeit  with  small  weighting  for  distance  samples)  from  the  entire  data  space  to  recon¬ 
struct  the  output  function  everywhere.  Of  course,  for  the  limited  data  record  such  as  in  SAU. 
the  global  interpolator  only  uses  data  from  the  limited  set  (outside  the  set.  samples  are 
assumed  to  be  zero.)  The  advantage  of  the  global  strategy  is  the  adherence  to  the  fact  that  a 
bandlimited  signal  is  analytic  (the  function  must  be  continuous  and  continuously 
differentiable)  and  thus  is  determined  everywhere  by  only  a  piece  of  the  function.  This  means 
that  when  interpolation  is  performed  in  the  transform  domain,  every  point  in  Fourier  space 
contributes  to  the  spatial  domain  reconstruction.  It  is  this  very  reason  that  low-order  interpo¬ 
lators  do  so  poorly  in  direct  Fourier  image  reconstruction  in  SAR  and  tomography. 

An  example  of  global  interpolation  is  the  zero  padding  of  a  DFT  sequence  prior  to  per¬ 
forming  the  inverse  FFT  to  obtain  a  finer  resolution  spatial  response.  This  is  often  used  in  con¬ 
junction  with  local  interpolators  to  interpolate  between  rotated,  rectangular  grids.  E.g..  to 
increase  image  resolution,  calculate  that  image  DFT.  zero  pad  by  a  power  of  2,  then  inverse 
DPT  the  zero  padded  sequence.  The  resulting  spatial  domain  resolution  will  be  increased  by  the 
padding  factor.  Nearest  neighbor  interpolation  (to  be  discussed  later  in  this  chapter)  is  then 
used  to  obtain  samples  on  the  rotated  rectangular  grid.  Speed  is  gained  by  means  of  the  FIT. 
though  at  the  expense  of  much  more  memory  usage.  In  a  2-1)  image  of  size  64  by  64  pixels 
(40%  data  points),  resolution  can  be  doubted  (in  both  x  and  y  directions),  but  memory  is 
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quadrupled  to  16384  data  points.  For  larger  data  arrays  (increased  resolution  factor),  the 
memory  required  may  be  prohibitive.  The  FlT-zero  pad  algorithm  necessitates  that  resolution 
be  uniformly  increased  everywhere  so  that  we  cannot  simply  increase  the  resolution  in  some 
small  neighborhood. 

It  is  important  to  place  the  zero  padding  in  the  correct  position  of  the  DFT  sequence.  The 
DFT  represents  uniform  samples  of  the  z-transform  of  the  sampled  signal  where  the  z- 
transform  of  the  finite  sequence.  h(n)  is  given  by 

H(z)  *  £  h(n)z"n  (3.10) 

n*<> 


j3wk 

Then,  substituting  z  -  e  N  .  we  obtain 


U(k)*  £  h(n)  e  ^  (3.11) 

««o 

The  dc  component  of  H(k)  is  located  at  k»()  (z  -  t)  and  the  highest  frequencies  are  at 
k-N/2  (z  -  -1).  The  zero  pad  therefore,  must  be  placed  at  both  ends  of  the  rotated  sequence, 
i.e..  in  the  middle  of  the  sequence  H(k).  This  is  at  the  high-frequency  folding  point  of  HCe^). 

Let  the  original  (time  domain)  sequence  be  x(n).  for  O^n^N— 1  (assume  N  even).  We 
wish  to  increase  resolution  by  a  factor  of  R  where  R  is  a  power  of  2.  Assume  also  that  N  is 
even,  since  most  applications  which  use  the  FFT  fix  N  to  be  a  power  of  2.  The  DFT  of  x(n)  is 

N-l  — 

X(k)*  r  x(n)  e  N  (3.12) 


The  dc  point  lies  at  X(0).  so  we  create  a  new  sequence  W(k)  by  splitting  X(k)  at  the  high- 
frequency  folding  point  and  inserting  I,  zeroes  where  1.  *  R  N: 


W(k)  = 


X(k) 

0 


X(L-k) 


0  <  k  ^  N/2  —  1 

N/2  k  <5  N/2  +  L  -  1 
N/2  +  L  <  k  ^  N+L-l 


(3.13) 


The  inverse  DFT  of  W(k)  will  yield  a  high  resolution  version  of  x(n) 


w(m)  =  N+I  W(k)  e+j^cri",k  (3.14) 

k=o 

This  FFT-zero  pad  algorithm  can  be  very  useful  for  determining  the  position  of  a  peak  in 
a  sampled  signal  when  it  falls  between  D1T  bins.  It  can  be  used  to  determine  the  movement  of 
a  point  target  from  one  range  (or  azimuth)  cell  to  another,  to  calculate  the  width  of  the  (possi¬ 
bly  moved)  peak,  or  to  align  sampled  waveforms  which  may  be  out  of  phase  [45].  As  an 
example,  consider  a  sampled  periodic  sine  waveform  (Figs.  3.2  and  3.3).  The  sine  peak  is  cen¬ 
tered  at  4.3.  but  the  sampled  sequence,  having  a  resolution  of  1.0.  could  be  misinterpreted  by 
assigning  the  peak  position  at  4.0.  The  resolution  can  be  doubled  by  calculating  the  DFT  of  the 
sequence  (F.q.  3.3),  zero  padding  the  sequence  to  2N.  then  inverse  transforming  the  result  (Fq. 
3.4).  After  this  first  resolution  enhancement,  the  peak  is  resolved  to  4.5.  This  can  be  bettered 
by  further  padding  the  DIT  sequence  to  4N.  8N.  and  16N.  each  time,  doubling  the  time  domain 
resolution  (Figs.  3.4  and  3.5).  In  one-dimension,  the  memory  usage  is  linear  with  the  interpola¬ 
tion  factor,  but  quadratic  in  2  D.  When  the  zero-pad  algorithm  is  used  to  interpolate  one  rec¬ 
tangular  grid  to  another  prior  to  nearest  neighbor  sampling,  the  memory  required  may  prohibit 
more  than  a  simple  doubling  of  the  sample  frequency. 

33  Separable  versus  Non-ae parable  Interpolators 


The  ideal  two-dimensional  reconstruction  kernel  is  a  bi-sinc,  or  a  separable  sine  function 
with  the  interpolation  kernel 


h(u.v) 


sin(~)  sin(^) 

*  U  *  V 


TTU 

Tu 


irv 

T7 


(3.15) 


This  is  easily  derived  from  the  one-dimensional  case  above,  and  like  the  one-dimensional  inter¬ 
polator.  suffers  from  the  same  problem  of  infinite  length  which  makes  it  difficult  to  implement. 
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More  important,  the  bi-sinc  is  only  separable  in  U  and  V.  not  in  R  and  0,  which  is  what  is 
required  for  the  polar  grid.  Thus,  what  is  normally  called  separable  for  the  Cartesian  grid  is 
really  non— separable  for  the  polar  grid.  The  separable  interpolator  is  one  that  can  be  decom¬ 
posed  into  two  stages.  The  first  stage  interpolates  the  polar  data  to  an  intermediate  grid,  from 
which  the  second  stage  interpolates  to  the  final  rectangular  raster.  The  first  stage  interpolates 
the  polar  data  to  a  keystone  grid,  and  the  second  stage  interpolates  from  the  keystone  grid  to  the 
final  rectangular  grid.  The  keystone  grid,  as  indicated  by  Fig.  3.6.  is  the  intersection  of  the  A0 
spaced  range  lines  and  AU  spaced  azimuth  lines.  When  the  range  lines  are  equi-spaced  in  0,  as 
in  a  standard  polar  grid,  the  azimuthal  spacing  on  the  keystone  grid  is  non-uniform..  If  the 
pulse  repetition  frequency  (PRF)  of  the  radar  is  constant,  then  the  azimuthal  spacing  is  con¬ 
stant  for  each  azimuthal  line  (lines  of  constant  U).  The  non-uniformity  of  the  azimuth  sam¬ 
ples  for  the  polar-keystone  grids  increase  the  amount  of  compulation  required  during  the 
second  stage  of  the  interpolation., 

3.4  A  Study  of  Interpolators 

It  is  clear  that  the  performance  of  an  interpolator  is  governed  by  its  approximation  to  the 
id?,*l  low-pass  filter.  It  is  important  to  have  very  high  rejection  in  the  stopband  and  be  as  fiat 
as  possible  in  the  passband.  Because  most  of  the  well-known  interpolators  are  of  limited 
(local)  support,  they  must,  by  definition,  be  suboptimal.  The  following  sections  describe  the 
spatial  domain  impulse  response  and  corresponding  Fourier  transform  of  several  well-known, 
but  rarely  characterized,  one-  and  two-dimensional  interpolators. 

3 .5  Nearest  Neighbor 

The  simplest  interpolator  is  the  zeroth  order,  or  nearest  neighbor.  It  produces  a  piecewise 
constant  function  (zeroth  order  polynomial)  and  is,  in  general,  discontinuous  at  each  break¬ 
point.  The  impulse  response  of  the  nearest  neighbor  interpolator  is  a  gate  (rectangle)  function 


Figure  3.6  Intermediate  Keystone  Grid. 


which  is  1  from  -T/2  lo  *T/2  end  0  outside  Urn  re*, on  (Fig.  3.7). 

h\'\(t)  a  pjy2(t) 

where 


fWO= 


1 

0 


III  <T/2 
III  >T/2 


(3.17) 


As  the  gate  function  is  convolved  with  thi»  innm  ria.  t 

tite  input  data,  exactly  one  input  datum  falls  within  the 

gate  and  will  be  selected  for  the  output..  Convolution 


in  the  spatial  domain  implies  multiplica- 


"I 
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Figure  3.7  Nearest  Neighbor  Impulse  Response. 


tion  of  the  transforms  in  the  frequency  domain.  The  transform  of  this  gate  function  is  a  sine 
waveform  with  very  high  sidelobes  (Fig.  3.8)  and  a  falloff  rate  of  3  db  per  octave. 


u  r-.'i-  2sin(a»T/2)  _  ,r  .  a»T 
Hnn(c«i)  = - - - *  T  sine  ~y 


(3.18) 


It  is  these  high  sidelobes  which  fail  to  remove  the  high-frequency  copies  of  the  original  func¬ 


tion.  It  is.  in  fact,  a  very  poor  low-pass  filter. 


The  nearest  neighbor  interpolator,;  as  defined  above,  is  also  ill-posed  for  points  lying 
exactly  half  way  between  known  data.  The  output  there  would  be  the  sum  of  the  two  neigh¬ 
bors.  This  can  be  averted  by  arbitrarily  chosing  only  one  of  the  points.  To  make  the  recon¬ 
structed  function  right-continuous,  we  chose  the  left  point  when  interpolating  exactly  at  the 
sample  midpoints. 
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Figure  3.8  Fourier  Transform  of  Nearest  Neighbor  Interpolator. 

In  two-dimensions,  the  nearest  neighbor  interpolator  is  a  2-D  gate-like  function  which  has 
a  shape  that  is  very  dependent  on  the  input  data  sampling  structure.  If  the  input  sample  spac¬ 
ing  is  rectangular,  then  the  nearest  neighbor  function  is  a  spatially  invariant  rectangular  gate 
pulse  (of  possibly  differing  u-v  dimensions).  If  the  input  array  is  not  uniformly  spaced,  then  it 
becomes  necessary  to  lessellate  the  input  sample  array  into  nearest  neighbor  regions  (Fig.  3.9). 
The  interpolator  then  becomes  spatially  varying  and  analysis  is  nearly  intractable  (in  the 
transform  domain)  except  possibly  in  a  stochastic  sense. 

On  a  regular  sample  grid  (e.g..  a  polar  grid)  it  is  possible  to  analyze  this  interpolator 
because  it  is  easier  to  parametize  the  2-D  gate.  For  a  polar  sample  grid,  the  nearest  neighbor 
function  is  a  spatially  varying,  rotated  trapezoid  which  grows  in  size  as  we  move  away  from 
the  origin,  and  which  rotates  with  the  9  sample  lines  (Fig.  3.10).  The  Fourier  transform  of  the 
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Figure  3.9  Tessellation  of  an  Irregularly  Spaced  Data  Array. 


rotated,  shifted  trapezoid  is  a  bit  messy,  but  straightforward.  The  spatially  varying  nature  of 
the  transform  results  in  a  transform  with  four  rather  than  two  parameters:  u.v  (frequency 
coordinates)  and  r.O  (spatial  coordinates)..  The  trapezoid  is  specified  by  three  parameters:  the  2 
bases  and  the  width  as  shown  in  Fig.  3.11.  The  Fourier  transform  of  this  trapezoid  (which  has 
the  value  1  inside  and  0  outside)  can  be  calculated  as  follows: 


First  we  see  that  the  trapezoid  is  really  a  linearly  warped,  shifted,  and  rotated  square.  The 
linear  warping  is  described  by  a  simple  transformation  matrix 


T(a  +  b)I  o 

-i(a2  +  b2)  £L(a+b) 


(3.19) 


i 
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This  matrix  transforms  a  square  centered  at  the  origin  of  width  one  to  the  trapezoid  centered  at 
the  origin  and  with  dimension  parameters  (a.b.c).  The  trapezoid  dimensions  (a.b.c)  are  depen¬ 
dent  on  the  spatial  position  of  the  trapezoid.  We  can  now  apply  a  linear  shift  transformation  in 
the  +u  direction  with  u’  -  u  +  A.  which  in  the  Fourier  domain  corresponds  to 

F(u.v)  =  e+,AuF(u,v)  (3.20) 

followed  by  a  rotation.  The  constant  A  represents  the  amount  of  spatial  translation  along  the 
U-axis.  A  is  a  function  of  the  position 


u" 

cosd  sind 

u* 

X 

v" 

—sind  cosd 

v' 

Hv  combining  (3.I9).(3.20).  and  (3.21).  the  spatially  varying  trapezoidal  impulse  response  and 
Fourier  transform  can  be  determined,  though  it  is  extremely  awkward  to  use  and  therefore  this 
formulation  has  only  nominal  value.  If  wc  convolve  this  spatially  varying  trapezoid  with  the 
input  polar  grid,  we  end  up  with  a  2-D  function  that  is  a  piecewise  constant  surface  that  looks 
like  blocks  of  varying  heights.  This  block-like  reconstruction  is  then  sampled  on  the  rectangu¬ 
lar  grid. 

A  problem  that  occurs  with  this  form  of  the  nearest  neighbor  interpolator  is  the  actual 
implementation.  Since  it  ,!s  not  easy  to  determine  which  trapezoid  the  rectangular  sample  falls 
on.  we  make  a  small  approximation  by  bending  this  shape  into  a  section  of  an  annulus  (Fig. 
3.12).  It  then  becomes  very  easy  to  determine  the  nearest  neighbor  with  simple  integer  round¬ 
ing.  The  index  of  the  angular  coordinate  is  found  by  adding  0.5  of  the  angular  sampling  spacing 
to  9.  and  then  truncating  to  the  integer  value.  The  range  index  is  found  in  a  like  manner  with 
the  range  sample  spacing.  For  the  particular  geometries  appropriate  in  spotlight  mode  SAR,  this 
is  a  very  close  approximation  to  the  ideal.  The  Fourier  transform  of  the  nearest  neighbor  inter¬ 
polator  is  shown  in  Fig.  3.13. 
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While  ihe  NN  algorithm  is  the  crudest  of  the  interpolators  presented  here,  it  still  requires 
square-root  and  arctangent  functions  to  locate  the  Cartesian  point  in  polar  space.  The  computa¬ 
tion,  however,  can  be  reduced  by  expanding  these  functions  in  a  Taylor  series  and  truncating  to 
3  terms.  The  approximation  is  good  enough  to  locate  the  point  in  polar  coordinates,  and  the 
error  in  position  is  usually  negligible  in  elTect,  i.e..  choosing  the  nearest  neighbor  incorrectly.  In 
our  simulations,  the  expansion  method  produced  incorrect  nearest  neighbors  in  only  3  out  of 
4096  points  -  0.073%  error.  Of  course,  the  square-root/arctangent  functions  could  also  be  gen¬ 
erated  through  a  lookup  table  combined  with  linear  interpolation.  As  we  shall  see.  these  small 
errors  in  the  nearest  neighbor  coordinate  approximations  are  usually  outweighed  by  the  error  in 
the  algorithm  itself. 

When  interpolating  from  one  rectangular  grid  to  another,  which  is  merely  displaced  in  U 
and  V  (no  relative  rotation),  where  the  input  and  output  rasters  have  the  same  sampling  fre¬ 
quency.  and  are  offset  from  one  another  by  AU.AV.  the  nearest  neighbor  will  merely  replicate 
the  original  data.  In  fad.  if  I  AC  I  <'I\V2  and  IAVI<Tv/2.  then  the  input  and  output 
responses  will  be  identical  [ll]. 

While  this  is  the  fastest  and  easiest  algorithm  which  can  be  used  to  produce  Cartesian 
samples  from  polar  formatted  data,  it  results  in  a  badly  distorted  response  for  CAT 
images  [24]  that  are  full  of  artifacts  and  false  lines.  In  SAR.  it  produces  many  false  targets. 

3.6  Linear  Interpolation 

Probably  the  most  widely  used  simple  interpolator  is  the  linear  interpolator.  It  is  com¬ 
monly  used  to  "read  between  the  lines"  of  tables  or  closely  spaced  samples.  When  a  function  is 
expanded  in  a  Taylor  series  about  some  known  point  £.  and  then  truncated  to  two  terms,  the 
resulting  expression  is  a  simple  linear  curve  (straight  line)  through  £.  A  signal  which  has  been 
reconstructed  with  a  linear  interpolator  will  be  continuous,  but  may  lack  continuous  first  and 
other  higher  order  derivatives  which  make  would  it  "smooth.” 
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The  linear  interpolator  is  usually  thought  of  as  a  connect-the-dol  interpolator.  It  is  more 
desirable  for  analytic  purposes,  however,  to  find  its  impulse  response.  This  is  the  convolution 
kernel  h(t)  in  our  original  reconstruction  formula 

v(t)  =  £  *(n)  h(  (3.22) 

|J  =  1 

where  T  is  the  input  sample  spacing.  One  formulation  of  the  linear  interpolator  is  given  by 
(3.23). 


*a-- 
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transforms,  resulting  in  a  sine  squared  Fourier  transform  (Fig.  3.14b). 

The  linear  interpolator  can  also  be  formulated  as  an  inverse  distance  interpolator  where 
y(t)  is  equal  to  the  weighted  sum  of  each  of  its  nearest  neighbors.  The  weights  are  the  normal¬ 
ized  reciprocals  of  their  distance  to  y(t),  i.e.. 


y(t)  = 


X(>+  x, 

l  T-t 

1+  « 

t  T-t 


(3.26a) 


=  x„(l  -  -1)  +  xi  ~ 


(3.26b) 


which  is  identical  to  (3.23b).  It  will  be  shown  that  this  relationship  does  not  hold  in  two- 
dimensions. 


3.6.1  One-dimensional  generalized  inverse  distance 

The  linear  or  inverse-distance  algorithm  can  be  generalized  to  the  inverse  distance  to  the  N 
(Nth  power  inverse  distance)  interpolator.  Here,  the  distances  to  each  of  the  two  nearest  data 
points  are  raised  to  some  power  N.  N>1.  This  results  in  weighting  the  nearer  point  more 
heavily  than  the  other.  It  also  results  in  a  continuous  first  derivative  in  the  reconstructed  signal 
at  each  interpolated  point.  The  impulse  response  of  the  generalized  inverse  distance  (G1D-N) 
interpolator  is 


h(t)  * 


(I  -d>s 
(1  -d)N  +  dN 


where 


(3.27) 


d  = 


III 

T 


This  is  obviously  not  a  polynomial  function,  but  rather  a  ratio  of  polynomials,  which  are 
more  difficult  to  analyze.  The  Fourier  transform  of  a  polynomial  is  easily  calculated,  but  the 
transform  of  this  class  of  functions  (defined  only  in  the  interval  [-T.T])  is  difficult  to 
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formulate.  However,  an  approximation  of  the  GID-N  transform  is  easily  accomplished  numeri¬ 
cally  with  the  DPT. 

In  the  ’imit  as  N  approaches  infinity.  GID-N  approaches  the  nearest  neighbor.  This  is 
obvious  because  the  weighting  of  the  nearer  point  approaches  infinity  (with  N)  and  causes  that 
point  to  be  selected  as  the  output.  With  very  higher  order  N.  T  must  be  scaled  to  prevent 
numeric  errors  due  to  computer  range  problems. 

3.7  Two-dimensional  Generalized  Inverse  Distance 


The  linear  interpolator  can  be  extended  to  two-dimensions  very  simply  for  rectangular 
input  grids.  For  this  case,  the  interpolator  can  be  thought  of  as  convolution  with  a  separable 
2-1)  triangle  signal  (a  quadratic  shaped  pyramid  (Fig.  3.15a)  shown  with  equal  spacing  in  u  and 
v)  having  a  half  width  equal  to  the  u  and  v  input  sample  spacings.  The  bilinear  interpolator 
has  a  Fourier  transform  which  is  a  separable  sine  squared  and  is  shown  in  Fig.  3.15b.  The 
impulse  response  for  the  bilinear  interpolator  is  given  by  (3.28). 


h(u.v)»(l  -  lul/Tu)(l  -  lvl/Tv)  lul  <Tu.lvl  <TV  (3.28) 


However,  when  the  input  grid  is  non-rectangular.  bilinear  interpolation  is  difficult  to  for¬ 
mulate  since  the  data  axes  are  not  orthogonal.  The  most  popular  linear  type  of  algorithm  used 
in  this  situation  is  the  inverse  distance  interpolator.  The  output  value  at  Q  is  the  weighted  sum 
of  the  4  nearest  neighbors  (Fig.  3.16). 


Q  = 


4  P 

TIL 

4  i 

rJL 


(3.29) 


The  impulse  response  is  given  by 


^  __  « 
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Figure  3.16  Inverse  Distance  Interpolator  Geometry. 


h(u.v)  = 


d2d3d4 

d2d3d4  +  d|d3d4  +  djd;jd4  +  djd2d3 


(3.31) 


1  ‘I 
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d4\ 

/d2 


h(u.v) 


ll/di 

»»i 


max(  I  u  I .  I  v  I )  <  1 


otherwise 


(3.30) 


where  d,  is  the  distance  from  (u.v)  to  (0.0).(0. ±  1 ).( *  1 ,0),  and  ( ±  1 . at  1 )  with  the  sign  depend¬ 
ing  on  the  quadrant  in  which  (u.v)  falls.  If  the  summation  in  (3.30)  is  expanded  and  the 
numerator  and  denominator  are  each  multiplied  by  d|d3d3d4,  h(u.v)  becomes 


which  avoids  the  problem  of  singularities  when  interpolating  at  an  input  data  point.  .Notice 
that  h(u.v)  has  the  value  1  at  (0,0)  and  becomes  0  at  the  eight  integer  sample  points  around  the 


origin  (0,0).(0,±  1 ),( ±  1 .0),( ±  1  ,±  1 ).  This  is  consistent  with  the  constraint  that  the  interpola¬ 
tor  should  exactly  reproduce  the  given  data  at  the  sample  coordinates.  The  inverse  distance 
impulse  response  and  corresponding  Fourier  transform  are  shown  in  Figs.  3.17a  and  3.17b. 

The  term  Generalized  Inverse  Distance  actually  stems  from  generalizing  the  order  of  the 
inverse  distance  weights  and  using  additional  local  points  in  the  summation.  These  points  are 
determined  by  extending  the  point  set  to  include  the  next  layer  of  points  around  the  nearest 
four..  The  first  additional  layer  contains  twelve  more  points  (sixteen  total),  the  second  layer 
contains  eighteen  points  beyond  that  (twenty-five  total).  This  admittedly  heuristic  point  set 
was  used  by  Shepard  [47]  in  contour  plotting.  It  is  used  here  in  the  computer  results  to  deter¬ 
mine  the  value  of  the  enlarged  point  set. 

It  was  demonstrated  in  the  II)  case  that  the  linear  and  inverse  distance  interpolators  were 
identical..  In  two-dimensions,  this  is  obviously  not  the  case,  since  (3.31)  cannot  be  reduced  to 
(3.28).  It  is  non-separable  and  difficult  to  analyze,  even  for  the  rectangular  grid  input  set.,  A 
comparison  of  the  transforms  for  the  bilinear  and  inverse  interpolators,  suggest  that  inverse 
distance  is  inferior  to  the  bilinear.  It  is  easy  to  see  (by  examining  the  first  derivative  behavior 
of  h(u.v)  about  the  (0.0)  point)  that  the  output  function  y(u.v)  will  not  have  continuous 
derivatives.  This  is  similar  to  the  bilinear  version,  though  more  pronounced,  due  to  the  steeper 
descent  from  (0.0).  The  problem  with  the  first  derivatives  can  be  rectified  by  using  inverse  dis¬ 
tance  squared  (G1D-N.  N-2).  This  choice  has  a  very  useful  computational  side  effect.: 

The  inverse  distance  squared  interpolator  has  some  very  nice  properties.  First,  the 
sidelobes  are  very  low  (Fig.  3.18b).  at  the  expense  of  a  slightly  wider  mainlobe. 


i=i 


For  the  SAR  application,  this  is  important  because  it  is  important  to  reduce  the  sidelobes  in 
order  to  reduce  the  artifacts  created  by  the  transform  operation.  Second,  the  inverse  distance 
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squared  algorithm  obviates  the  need  to  perform  the  square-root  operation  in  calculating  the  dis¬ 
tances  to  the  nearest  4  points.  This  is  quite  a  computational  savings,  as  will  be  shown  in  the 
computer  simulalicas.  Surprisingly,  this  algorithm  has  had  little  application  in  the  image  pro¬ 
cessing  field,  but  rather  has  been  used  more  in  such  areas  as  contour  and  2-D  surface  plotting 
(graphics)  [47],  Notice  that,  like  the  one-dimensional  case,  as  N  approaches  infinity,  the  G1D-N 
interpolator  approaches  the  NN  interpolator  (Figs.  3.19  and  3.20). 

3.7.1  Arbitrary  local  point  weighting 

In  the  most  general  case,  when  interpolating  to  a  point  Q.  the  nearest  1*  points.  ps.  are  indi¬ 
vidually  weighted  by  a  function  w,  a  function  of  the  Cartesian  distance  from  Q  to  P,„ 
d,  [47.  48].  That  is. 

IW(d.)Tl 

Q  =  — -  (3.33) 

£w(dt) 

IN 

The  set  of  points  p,  is  usually  chosen  to  be  the  nearest  P  points,  or  a  set  of  point  within  a 
given  circle  (used  in  irregularly  space  data  sets  so  that  areas  of  sparse  data  do  not  used  very 
distance  elements  in  the  weighted  sum.)  w(d)  can  be  any  monotonically  decreasing  function 
which  satisfies  (3.34a)  and  (3.34b). 

lim  w(d)  =  0  (3  34a) 

d-*oo 

lim  w(d)  =  oo  (3.34b) 

This  very  simple  constraint  will  reproduce  the  sampled  data  exactly  and  can  lead  to  several 
heuristic  weighting  functions  such  as 

w(d)=— i_  (3.35) 

e" — 1 

which  has  an  exponential  weighting  surface.  In  the  previous  section.;  w(d)  =  d".,  where  n  =  1 
for  inverse  distance  weighting  and  n  =  2  for  inverse  distance  squared. 
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3.7.2  Warped  and  spatially  varying  nature  of  GED 

It  should  be  pointed  out  that  the  Generalized  Inverse  Distance  interpolator  suffers  from 
the  same  spatially  varying  problem  that  the  nearest  neighbor  algorithm  does.  For  a  polar  grid, 
the  base  of  the  tent  function  (for  simple  inverse  distance)  grows  larger  and  more  warped.  In 
the  Nearest  Neighbor  section,  the  frequency  response  was  generated  by  warping  (and  rotating) 
a  square  gate  pulse.  This  could  be  done  for  the  inverse  distance,  but  little  insight  would  be 
gained  by  such  mathematical  manipulations.  Rather,  the  very  narrow  look  angle  typical  of 
spotlight  mode  SAR  allows  us  to  approximate  the  warped,  spatially  varying  impulse  response 
with  those  given  above  (shown  for  square  input  spacing)  as  a  stationary  function. 

3.8  The  Weighted  Sine  Interpolator 

A  third  interpolator,  henceforth  to  be  called  the  Weighted  Sine  interpolator  (also  referred 
to  in  the  literature  as  the  Windowed-Sine.  Standard  Polar  Format,  and  ERIM  interpolator),  w- 
sinc.  is  a  two-step  algorithm.  First,  the  polar  samples  are  interpolated  along  the  range  lines  to  a 
"keystone  raster."  that  is.  a  sampling  raster  which  has  data  on  parallel  azimuthal  lines  (Fig. 
3.21).  Note  that  this  interpolation  is  implemented  using  one-dimensional  algorithms  on  a  range 
line  by  range  line  basis.  In  a  sense,  it  is  separable,  though  not  in  the  typical  U-V  coordinate 
system,  but  rather  in  the  R.0  coordinate  space.  In  step  two.  the  data  are  interpolated,  again 
using  one-dimensional  techniques,  to  the  final  rectangular  grid.  The  computational  complexity 
of  this  algorithm  is  determined  by  the  complexity  of  the  1-D  interpolator  for  each  step. 

It  was  shown  earlier  that  the  sine  function  is  the  ideal  reconstruction  filter,  although  it  is 
impractical  due  to  its  infinte  support.  If  we  attempt  to  use  a  truncated  version  of  the  sine,  the 
resulting  transform  displays  the  usual  Gibb’s  phenomenon  oscillations  which  result  from  high 
sidelobes  and  poor  passband  response.  This  problem  can  be  alleviated  by  truncating  the  sine 
function  with  a  tapered  window  function. 


Figure  3.21  Interpolation  to  an  Intermediate  Keystone  Raster. 

sin(  ) 

h(t)  *  w(t)  • - (3.36) 

w  l 
"Tout 

where  rou,  is  the  sampling  interval  ol  the  output  sequence  and  w(t)  is  the  windowing  function., 
I  he  sine  1  unction  is  specified  in  terms  of  the  output  sample  "pacing  because  it  is  actually  doing 
the  job  of  low-pass  filtering  and  signal  reconstruction  in  one  convolution.  The  cutoff  frequency 
of  the  low-pass  filter  is  given  by  the  output  spacing  so  that  when  we  interpolate  to  a  sparser 
array  ( F>T),  the  filtering  operation  will  prevent  aliasing  of  high-frequency  data  into  the  low 
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frequencies.  In  the  first  stage  of  interpolation  from  the  polar  grid  to  a  keystone  grid,  the  input 
sample  spacing,  AR,  is  constant  from  range  line  to  range  line,  but  the  output  sample  spacing 
changes,  i.e..  the  sine  width  changes  very  slightly  for  each  range  line  interpolation. 

Windows  are  most  typically  used  to  taper  data  records  prior  to  applying  a  transform 
(either  DFT  or  IDFT)  to  reduce  the  sidelobes  caused  by  truncating  the  infinite  sample  sequence. 
Harris  [49]  has  cataloged  a  large  list  of  well-known  windows  and  their  spectrum  shaping  pro¬ 
perties  in  conjunction  with  spectral  estimation  and  spectral  leakage.  Although  many  windows 
are  candidates  for  the  w-sinc  interpolator,  the  Hamming  window  was  chosen  in  the  computer 
simulations  because  of  the  ease  of  evaluation  and  familiarity.  The  Hamming  w-sinc  interpola¬ 
tor  becomes 

sin(  —■ ) 

h(t)=  0.54  +  0.46  cos(i£!)  • - —  (3.37) 

K  n  t 

or  more  simply 

h(t)=  0.54  +  0.46  cos(  )  sincf-^)  (3.38) 

where  K  specifies  the  support  of  the  sine  calculated  in  terms  of  the  output  sample  spacing. 
Because  of  this  specification  of  K.  the  number  of  input  samples  that  fall  within  h(-)  may  vary 
by  1  between  any  two  interpolator  sums.  This  has  negligible  effects. 

The  Hamming  window  is  actually  a  modified  Hanning  window  which  is  a  member  of  the 
cos“(0)  class  (a  *  2  for  Hanning).  It  can  be  formulated  as  the  multiplication  of  a  raised  cosine 
by  a  rectangular  window  -  in  the  transform  domain,  it  is  the  convolution  of  the  sine  kernel 
with  three  impulses  located  at  0  (of  height  2ir).  and  at  ±  it/ K  of  heights  ir  each.  Its  transform, 
therefore.  Is  the  weighted  sum  of  three  sine  kernels  [49] 

W (<u)  =  0.54'2rr  sinc(^-w)  +  0.46-ir  sincf-JL  +  o>))  +  0.46-7T  sinc(^-  —  tu))  (3.39) 
K  K  K 

The  side-lobes  of  the  two  offset  sines  cancel  out  the  sidelobes  of  the  main  sine,  and  thus,  the 
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overall  transform  has  lower  sidelobes  at  the  expense  of  the  wider  main  lobe.. 

The  Fourier  transform  of  the  Hamming  windowed  sine  is  given  as  the  FT  of  h(l): 

H(w)  =  FT  I  sinc( -Z-  t)  -  (0.46  +  0.54  cos(  Z  l))  *  pK(t)  }  (3.40) 

*  out  k 


which  is  the  convolution  of  the  transforms  of  each  term.  The  transform  of  the  Hamming  win 
dow  is  given  in  (3.41). 


(I.OXir2  —  0.16K2<t>2)sin(TTto) 
aAir2  —  K2o>2) 


(3.41) 


As  we  increase  K.  we  take  more  terms  into  the  summation.  This  has  the  effect  of  creating 
a  better  low-pass  filter,  i.e..  better  interpolator,  but  also  of  increasing  processing  lime  linearly 
with  K.  Figures  3.22  and  3.23  show  the  impulse  response  and  corresponding  Fourier  transform 
of  the  w-sinc  for  K-4  and  k-10.  While  it  is  obvious  that  the  w-sinc  shows  the  most  promise 
for  an  interpolation  kernel,  the  drawback  is  the  excessive  amount  of  computation  required  in 
the  lengthy  summation  (direct  convolution),  the  evaluation  of  the  sine  and  cosine  functions, 
and  the  calculation  of  the  input  data  positions.  The  data  positions  of  an  a  priori  known  polar 
grid  could,  in  theory,  be  stored  in  a  ROM:  however,  the  data  collection  path  of  a  maneuvering 
aircraft  precludes  such  a  ROM.  Hence,  data  positions  must  be  calculated  "on  the  fly  "  Appendix 
11  describes  a  novel  method  of  reducing  sinusoid  compulation,  but  the  w-sinc  interpolator  is 
still  costly. 


3.9  Splines 


3.9.1  Polynomial  interpolation 

The  use  of  polynomials  for  function  approximation  and  interpolation  is  widespread 
because  polynomials  are  very  easy  to  manipulate,  differentiate,  and  integrate.  They  are  ana¬ 
lytic  and  well-behaved  and  are  representable  in  a  large  number  of  ways.  The  Taylor  expansion 
of  a  function  is  given  as  a  (possibly  infinite)  polynomial  which  is  often  approximated  by  trun- 


Figure  3.22a  Impulse  Response  of  Windowed  Sine  of  Order  4 

FREQUENCY  RESPONSE 


Figure  3.23a  Impulse  Response  of  Windowed  Sine  of  Order  10. 

FREQUENCY  RESPONSE 


NORMRUZEO  FREQUENCY 

Figure  3.23b  Fourier  Transform  of  Windowed  Sine  of  Order  10, 
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eating  to  a  finite  length  polynomial.  In  computer  applications,  polynomials  are  easily 
represented  as  an  array  of  coefficients  (or  derivatives  at  a  specific  point).  The  linear  interpola¬ 
tor  presented  in  a  prior  section  is  really  a  first  degree  polynomial6.  It  is  not  surprising,  there¬ 
fore,  that  polynomials  could  be  used  for  our  frequency  domain  interpolation  problem. 

Given  a  set  of  data  points  yj.yj. .y^  at  lime  samples  Xj.Xt.  • •  .Xyj.  it  is  well  known 
that  an  N-l  th  degree  polynomial,  pyj(x).  can  be  generated  such  that  Ps(Xj)  =  yp  T 'he  Nth  order 
polynomial  can  be  written  as: 

p(x)  =  at  +  a2x  +  •  •  •  +  a„x''-1  (3.42a) 

=  £a,xi-'  (3.42b) 

I«1 


where  the  coefficients  a,  can  lie  generated  by  first  forming  the  intermediate  polynomials  (called 
the  LaGrange  polynomials) 


N  v  — . 


il(X)= n-^ — 

7»»x*“xj 


(3.43) 


which  have  the  property  that 


l,(Xj)  = 


0  i^j 
1  i=j 


(3.44) 


The  nth  order  polynomial  is  then  generated  through  the  sum  of  these  LaGrange  polynomials 

p(x)  =  £  >T  (3.45) 

i=i 

Note,  however,  that  this  becomes  rather  unwieldy  to  work  with  for  large  N,  It  is  also  computa¬ 
tionally  unstable  on  a  finite  precision  computer. 


6 Some  texts  call  it  4  second  order  polynomial  because  it  is  specified  by  two  coefficients  -  it  spans  4  two-dimensional 
space.  The  terms  degree  and  order  shall  be  treated  here  as  the  highest  power  ol  the  polynomial,  i.e.,  a  cubic  is  a  third 
degree  polynomial. 


Another  problem  with  polynomials  of  increasing  degree  is  convergence.  As  the  number  of 
interpolating  points  is  increased,  the  interpolating  polynomial  pu(x)  is  not  guaranteed  to  con¬ 
verge  to  the  original  function  f(x).  In  fact,  it  can  be  shown  that  for  any  sequence  of  sets  of 
points  x,.  there  is  a  function  f(x).  such  that  the  sequence  of  interpolating  polynomials  will 
diverge  at  these  points  [34].  Finally,  the  multidimensional  polynomial  is  not  unique  for  a 
given  set  of  points,  i.e..  a  rectangular  grid  of  sample  points  will  not  produce  a  unique  polyno¬ 
mial  of  the  form 

p(x.y)  =  a„  x"  +  a„_|  x"-1y  +  a„_2  x'1-2y2  +  *  •  •  +  a,  y"  +  a»  .  (3.46) 

This  means  that  it  cannot  even  be  used  as  a  local  2-D  interpolator,  precluding  its  use  for  the 
SAR  system. 

Though  there  are  various  ways  to  represent  a  polynomial  which  can  reduce  the  amount  of 
instability  in  computation  (Newton’s  form  with  divided  differences,  for  instance),  these  Nth 
order  functions  tend  to  oscillate  too  much  between  the  nodes.  This,  of  course,  is  because  of  the 
fact  that  an  Nth  order  polynomial  must  have  N  roots.  On  a  digital  computer,  the  high-order 
terms  of  the  polynomial  are  very  sensitive  to  small  changes  in  x.  This  is  sometimes  referred  to 
as  an  ill-conditioned  function.  One  well  known  technique  to  deal  with  these  shortcomings  is  by 
using  different  lower  order  polynomials  in  different  regions  of  the  curve.  The  result  is  the  class 
of  polynomials  called  piecewise  polynomials  or  pp  functions.  Another  name  for  pp  functions  is 
splines , 

The  term  spline  has  a  rather  broad  definition  that  is  frequently  misused  in  the  literature. 
Most  simply,  a  spline  is  a  piecewise  polynomial  without  constraints,  i.e..  nothing  is  said  about 
continuity  of  the  spline  or  continuity  of  any  of  its  derivatives.  Typically,  however,  the  spline 
is  constrained  to  be  at  least  continuous,  and  m  derivatives  of  the  Nth  order  spline  are  to  be 
continuous,  where  m  <  N— 1.  The  linear  interpolator  can  be  referred  to  as  a  first-degree  spline 
which  is  continuous,  but  lacks  any  continuous  derivatives.  A  parabolic  spline  is  a  piecewise 
quadratic  function  which  has  a  continuous  first  derivative.  A  cubic  spline  is  continuous  and  has 
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continuous  first  and  second  derivatives  everywhere  on  the  interval. 

Spline  functions  are  often  used  in  mathematical  and  curve-fitting  applications,  because 
they  can  generate  smooth  interpolating  functions  with  continuous  derivatives  at  the  known 
data  points  (also  called  “knots."  “joints.”  and  "breakpoints")  [50,51,52].  More  recently, 
splines  have  been  used  in  computer  graphics  to  generate  smooth  curves  from  discrete  edge 
points  of  objects,  and  thus  reduce  the  number  of  points  required  for  an  irregular  curve. 

Splines  have  also  been  used  in  typesetting  applications  where  enlargement  and  reduction 
of  character  sets  are  more  easily  accomplished  with  pen  (plotter,  laser,  light)  strokes  than  with 
binary  digitizations  [53].  By  storing  character  sets  as  a  combination  of  straight  lines  and  spline 
curves,  they  can  be  arbitrarily  scaled  without  degradation, 

3.9.2  Cubic  spline  interpolation 

To  see  how  splines  can  be  constructed  to  be  smooth,  i.e..  to  possess  first  and  second  deriva¬ 
tives  at  each  knot,  we  construct  a  simple  case  with  only  two  knots  with  known  values  and 
derivatives.  The  following  is  taken  from  Rivlin  [34]. 


If  ot<(i  .  and 

I 

p(a)  =  Uj  p(0)  -  u2 
p'(o)  *  u'j  p'(0)  -  u\ 


(3.47) 


then 


p(x)  =  u, 


(x  -  0)2 

(0  -  a)2 


j_  (x  —  a)(x  —  0)2 
(0  —  a)3 


(x  —  a)2  _  2(x  —  0)(x  —  a)2 
(0  —  a)2  (0-a)3 


X  ...  (x-a)(x-0)2  ^  (x  -  >x)2(x  -  0) 
(0  —  a)2  *  (0  —  a)- 


(3.48) 


Because  a  cubic  polynomial  has  4  degrees  ol  freedom,  the  polynomial  can  be  forced  to  have  a 
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given  value  and  derivative  at  each  of  2  points. 

If  the  spline  is  constrained  to  have  a  continuous  second  derivative,  then  it  can  be  shown 
that  the  entire  spline,  s(x),  which  consists  of  the  cubic  polynomial  sections, 
s,(x)  for  x,^x<x1+1,  can  Ire  uniquely  specified  by  the  given  breakpoint.;  Xj,  the  sample  values  y,, 
and  the  boundary  conditions  u'j  and  u'N.  The  system  of  equations  thus  depends  on  the  initial 
and  final  derivative  of  the  knot  sequence.  The  end  conditions  can  also  be  specified  as  second 
derivati  .  >\  or  constraints  on  the  first  and/or  second  derivatives  which  simply  increase  the  sys¬ 
tem  of  equations  which  must  be  solved  (by  2  more  equations).  It  is  interesting  to  note  that 
spline  interpolation  has  a  global  nature.  Each  datum  in  the  set  of  N  points  affects  the  entire 
spline  function,  though  in  a  much  more  subtle  way  than  the  Nth  order  polynomial  interpolant. 
Delloor  showed  that  this  system  of  equations  forms  a  band  symmetric  matrix  which  can  be 
inverted  by  Gaussian  elimination  without  pivoting  [54].  The  recursive  solution  is  very  fast  and 
is  implemented  in  the  IMSL  mathematics  programming  library. 

The  various  boundary  conditions  which  can  be  imposed  on  the  spline  function  are  dis¬ 
cussed  here.  Unfortunately,  the  names  of  the  more  popular  ones  are  not  wholly  descriptive  of 
their  effects.  If  the  derivative  of  the  original  function  is  known  at  the  end  points,  then  we 
select  u'jfx,)  *  f*(x,)  and  u'n(x,)  *  f(xN).  The  resultant  spline  is  called  the  complete  cubic 
spline  of  f.  Most  DSP  applications  do  not  have  derivative  values  available,  nor  are  they 
estimated  well  due  to  additive  system  noise.  If  the  second  derivatives  at  the  endpoints  are 
forced  to  be  0.  then  s(x)  becomes  the  natural  spline.  This  is  a  popular  boundary  condition  in 
available  software,  but  is  rather  arbitrary  and  leads  to  only  (X  III2)  convergence,  which  means 
that  as  we  decrease  the  spacing  between  the  knots,  the  function  s(x)  converges  to  f(x)  in  a 
squared  fashion.  By  contrast,  the  complete  spline  interpolant  converges  in  0(  It!4)  to  f(x). 

Another  possible  boundary  condition  is  called  not— a— knot.  In  this  case,  the  first  and  last 
knots  of  the  sequence.  xt  and  xN.  become  inactive  and  sj(x)  =  s2(x)  and  sN_j(x)  =  sN_2(x).  This, 
too.  has  (X  1 1 1  4)  convergence  to  f(x).  It  is  the  end  condition  used  in  the  I  MSI,  math  library  for 
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spline  interpolation. 

The  cubic  spline  function  is  given  by  piecewise  cubic  polynomials  of  the  following 
form  [54].: 

s,(x)  =  c, ,  +  c2.,  (x  -  I,)  +  c3l  (x  - 1,)2  +  c4l  (x  -  D3  x,  <  |,  <  x1+,  (3.49) 

where  the  coefficients  c,  j  are  chosen  such  that 

s,(x,)  =  s,_j(x,)  l^i^N  (3.50a) 

s',(x,)  =  s'j _ ,( x,)  2<  i  <N-1  (3.50b) 


The  constant  coefficient  is  set  to  ct ,  =  v,  so  that  the  spline  interpolates  to  the  input  data  set  at 
the  knot  positions.  Once  these  coefficients  have  been  determined  through  the  Gaussian  elimina¬ 
tion  step,  it  is  a  simple  matter  to  locate  the  interpolation  point  within  the  correct  interval  and 
compute  the  cubic  polynomial.:  If  many  points  are  to  be  generated  within  an  interval,  this 
method  would  be  more  preferable  than  the  weighted  sine  from  a  computational  view. 

The  cubic  spline  interpolant  has  the  unique  property  of  maximal  smoothness.  That  is,  for 
a  given  set  of  interpolation  points,  the  cubic  spline  gives  the  smoothest  curve  over  all  functions. 
The  smoothness  is  defined  as  the  integral  of  the  square  of  the  second  derivative: 

t> 

/[f"(t)]2dt  (3.51) 


This  integral  is  minimized  over  all  f(t)  by  s(t).  the  spline  interpolant.  It  is  also  termed  the 
roughness  [55].  The  smoothness  integral  is  an  approximation  to  the  strain  energy  function 


f  [r(t)]2dt 

{  [1  +  (f'(l))2F2 


(3.52) 


This  is.  in  fact,  where  the  term  spline  is  derived.  A  spline  was  a  thin  rod  used  by  draftsman  to 
fair  curves  through  a  set  of  data  points.  The  curve  naturally  minimized  its  strain  energy  and  is 
approximated  by  the  cubic  spline  function. 
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Because  the  spline  interpolating  function  is  global  and  requires  a  matrix  solution,  it  does 
not  possess  a  time  invariant  impulse  response.  Horowitz  has  determined  the  power  spectral 
density  effects  of  spline  reconstruction  [56].  He  found  that  higher  order  splines  are  needed  to 
preserve  the  power  spectrum  of  the  original  function,  but  that  these  higher  order  polynomials 
require  much  more  compulation  to  generate  and  evaluate.  A  close  cousin  to  the  complete  spline 
interpolant  is  formed  via  the  B-spiine  functions. 

3.9.3  B-splines 

Hou  and  Andrews  [57]  used  cubic  basis  splines,  or  B-splines  (sometimes  called  0-splines). 
for  image  enlargement  and  reduction,  and  for  text  magnification  and  minificalion.  with  the  con¬ 
clusion  that  the  B-spline  may  prove  more  appropriate  than  a  truncated  sine  for  finite  length 
data  records. 

The  term  B— spline  comes  from  the  basis  formulation  of  the  reconstruction  function.  If 
f(t)  is  the  original  sampled  function  (to  be  reconstructed)  at  the  points  x(nT).  then  we  con¬ 
struct  a  set  of  basis  functions,  B,(t),  and  form  the  sum 

f(t)  =  £c,B,(l)  (3.53) 

i=i 

Note  that  the  c,  are  not  the  original  samples  y,.  but  rather  new  coefficients  calculated  from  y,. 
For  equally  spaced  data,  the  B-spline  functions  are  shift  invariant  and  the  subscript  i  can  be 
dropped  f rom  B. 

The  cubic  B-spline  is  formally  defined  as  the  convolution  of  two  triangle  (Chateau  func¬ 
tions)  pulses.  ea<;h  of  width  T  and  heigl  t  1,  The  resulting  function  is  given  as  two  cubic  poly¬ 
nomials  and  is  symmetric  about  the  origin.  The  interpolation  keinel  is  given  by:: 

l3/2  - 12  +  2/3  i  1 1  <  l 

h(t)=  — 13/6  +  t2  -  2t  +  4/3  l<ltl<2  (3.54) 

0  Itl  >  2 


X. 


w 
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This  can  be  recognized  as  the  Parzen  window  which  has  a  Fourier  transform  (scaled  so  that 
F(0)=  1). 


F(ta)  = 


sinc( 


oi  T 
~T 


) 


i* 


(3.55) 


This  is  expected  since  the  convolution  of  two  triangle  waveforms  corresponds  to  the  product  of 
its  transform,  which  are  in  turn  resulting  from  the  convolution  of  two  gate  pulses.  The  nth 
order  B-spline  can  be  recursively  generated  by  convolving  the  (n-l)th  order  B-spline  with  a 
gate  pulse  of  width  T. 

This  formulation  of  spline  reconstruction  is  actually  a  different  interpretation  of  the  com¬ 
plete  spline  given  above.  The  implementation  is  slightly  different,  but  the  dependency  of  c,  on 
y,  is  global.,  and  thus  while  the  B-spline  does  have  finite  support,  the  B-spline  coefficients  do 
not.  The  interpolation  function  is  generated  by  computing  the  c,'s  and  then  filtering  this 
sequence  with  the  B-spline  kernel.  The  previous  interpretation  leads  to  a  faster  implementation 
and  is  used  in  the  computer  analysis.  A  B-spline  type  of  reconstruction  can  be  done  with  a 
slight  modification  to  the  basis  function,  which  leads  to  an  interpolation  impulse  response  that 
resembles  a  truncated  sine  function.. 


3.9.4  Cubic  spline  convolution 

The  complete  spline  interpolant  is  generated  by  solving  a  system  of  equations  which  come 
about  by  the  constraints  on  the  interpolating  points,  the  first  and  second  derivative  matching  at 
these  points,  and  the  boundary  conditions.  A  more  general  cubic  spline  that  does  not  have  all 
these  constraints  is  given  by  the  following  piecewise  cubic  polynomial  which  is  symmetric 
about  the  origin  [  1 1  ]: 


a3t3  +  a2t2  +  ajt  +  a0 
b3t3  +  b2ti  +  bjt  +  b„ 


h(t) 


Itl  <1 
1  <  Itl  <2 


(3.56) 
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To  create  an  impulse  response  that  interpolates  the  data  exactly.  f(x)  must  be  0  at  t=l  and 
t=2.  Also,  first  derivative  of  each  cubic  should  match  at  1=1..  This  gives  seven  constraints  in 
all.  but  the  piecewise  function  has  eight  constants.  This  constant  can  be  made  a  variable  in  the 
function  f(x)  which  can  be  modified  to  generate  different  cubic  spline  functions.  It  is  important 
to  remember  that  this  spline  is  different  from  the  complete  cubic  spline  interpolant  presented 
previously. 


h(t) 


(c+2)t3  —  (c+3)t2  +  1  1 1 1  <  1 

cts  — •  5ct2  +  Set  —  4c  1  ^  Itl  ^2 


(3.57) 


When  the  constant  c  has  the  value  -0.5.  the  interpolation  error  goes  to  zero  with  the  third 
power  of  the  sample  interval  (OIT3l  ).  The  second  derivatives  of  the  two  polynomials  are 
forced  to  be  equal  if  c  =  -0.75.  This  same  result  was  used  by  Keys  [58]  who  demonstrated  that 
the  accuracy  of  the  interpolation  is  between  the  linear  and  complete  spline  interpolant.  This 
form  of  cubic  spline  convolution  is  studied  in  the  computer  analysis  of  Chapter  5.  The  impulse 
response  of  the  modified  B-spline  with  c  -  -0.5  is  plotted  against  the  w-sinc  of  order  4  in  Fig. 
3.24a.  Note  that  this  B-spline  is  remarkably  similar  to  the  w-sinc.  Similar  performance  is 
expected.  The  Fourier  transform  of  the  B-spline  (c  •  -0.5)  is  presented  in  Fig.  3.24b. 


COMPARISON  OF  S1NC  (4) 


BSPLINE  (-0.5) 


ANO 


Figure  3.24a  B-spline  (c  «  -0.5)  Impulse  Response  with  Sine  (order  -  4). 

FREQUENCY  RESPONSE 


Figure  3.24b  Fourier  Transform  of  B-spline  (c 
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CHAPTER  4 

WINDOWS,  SPURIOUS  TARGETS  AND  ALTERNA  FIVE  SAMPLING  GRIDS 

This  chapter  presents  a  number  of  topics  that  are  important  in  SAR  processing.  Window 
design  for  SAR  data  arrays  is  discussed,  and  the  usefulness  of  optimal  windows  is  examined.. 
Spurious  target  generation  resulting  from  the  Pourier  domain  interpolation  stage  is  presented, 
and  some  one-dimensional  examples  are  shown.  Finally,  alternate  sampling  grids  are  proposed 
which  can  reduce  the  computation  and  even  the  complexity  of  the  sampling  hardware. 

4.1  Windowing  the  Data  Array 

It  has  been  demonstrated  that  spotlight  mode  SAR  can  generate  a  data  set  which  is 
approximately  the  Fourier  transform  of  the  complex  reflectivity  map.  This  is  inverted  with  the 
DFT  (FFT).  The  same  procedure  is  used  for  the  tomographic  reconstruction  problem  as  well  as 
remote  sensing  arrays  (ground-based  aeronomy  SARs).  However,  when  using  the  DFT  as  an 
approximation  to  the  Fourier  integral,  errors  occur  due  to  the  limited  .aze  and  particular  shape 
of  the  recorded  data.  It  is  well-known  that  strict  truncation  (uniform  windowing)  of  a  signal 
results  in  Gibb’s  oscillations,  or  sidelobes.  which  can  obscure  weak  signals.  These  sid  Jobes  can 
be  reduced  through  the  use  of  a  weighting  function,  or  window  which  tapers  to  zero  near  the 
ends  of  the  sampled  data  set.  In  addition,  this  tapering  smooths  out  discontinuities  at  the  end 
points  (the  DFT  assumes  that  the  data  represents  one  period  of  a  periodic  signal,  so  x(O)  follows 
x(N-l))  which  can  result  in  sidelobes  as  well.  The  type  of  window  used  dramatically  affects 
the  resulting  output  spectrum,  often  reducing  the  high  sidelobes  at  the  expense  ol  a  wider 
mainlobe.  The  disadvantages  of  windows  is  the  loss  of  inlormation  from  the  tapered  spectral 
shaping  which  may  not  be  representative  of  the  signal  spectrum,  and  a  more  difficult  evaluation 
of  the  output  image,  i.e..  the  MSI:  cannot  be  used  as  a  good  measure  of  image  quality.  Win¬ 
dows  have  also  been  used  success)  ully  i  the  design  of  FIR  filters  [59].  Indeed,  the  weighted 
sine  interpolator  ol  Chapter  3  was  essentially  a  filler  designed  by  tapering  the  sine  function 
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with  a  Hamming  window. 

Harris  [49]  has  examined  a  number  of  windows  and  cataloged  them  according  to  peak 
sidelobe  levels,  mainlobe  width,  energy  leakage,  and  several  other  parameters.  The  windows  are 
also  presented  graphically,  and  examples  of  using  windows  in  spectral  estimation  are  given. 
Additional  window  designs  and  measurements  are  given  by  Geckinli  and  Yavuz  [60]  which 
favor  the  Hamming  and  Kaiser-type  windows  for  spectral  analysis  and  demonstrate  the 
optimality  of  the  truncated  sine  window  for  an  energy  maximization  criterion.  Because  SAR 
can  be  viewed  as  a  spectral  estimation  problem,  windows  are  examined  which  can  improve 
sidelobe  reduction  and  resolve  point  targets. 

In  this  section,  we  briefly  examine  some  of  the  windows  which  were  used  in  the  computer 
evaluations  of  2-1)  interpolators.  Previous  work  in  optimal  window  design  for  irregularly 
shaped  data  set  is  also  discussed. 

4.1.1  Separable  versus  circularly  symmetric  windows 

Probably  the  most  common  windows  used  in  two-dimensional  signal  processing  are  separ¬ 
able.  That  is.  they  can  be  represented  as 

w(x.y)  *  w,(x)  •  w2(y) 

where  typically.  w(  and  w3  are  the  same  function,  which  is  a  good,  one-dimensional  window. 

Huang  [61]  has  shown  that  good  circularly  symmetric  filters  can  be  designed  by  rotating 
good  one-dimensional  windows,  i.e.. 

w2_D(x.y)  =  wu,(V x2  4-  y2) 

The  circular  weighting  function  has  been  using  in  optical  SAR  processing  due  to  the  nature  of 
the  spherical  lens.  It  was  also  the  weighting  function  chosen  by  ERJM  for  their  processing  of  a 
rotating  platform  imaging  system  [62]. 


Unfortunately,  when  using  windows  with  the  DPT  for  spectral  analysis,  a  window  which 
is  non-zero  (although  equal)  at  the  endpoints  will  generate  discontinuities  when  rotated,  and 
thus  produce  higher  sidelobes  than  the  one-dimensional  case.  A  Hamming  window  is  one  exam¬ 
ple.  The  two-dimensional  circularly  symmetric  Hamming  window  shown  in  Fig.  4.1b  displays 
the  discontinuous  boundary.  This  is  an  example  of  a  window  which  is  good  for  the  one¬ 
dimensional  and  separable  2-1)  cases,  but  possesses  a  poor  circular  response.  Experimental 
results  in  the  next  chapter  demonstrate  the  inferior  image  quality.  This  may  be  corrected  with 
a  zero  taper  window  such  as  the  Hanning  window  which  displays  no  such  discontinuities  (Fig. 
4.2).  The  additional  problem  with  the  rotated  window  is  the  great  information  loss  due  to  the 
relatively  large  area  of  zero  weighting  near  the  high  frequency  "corners."  If  the  energy  concen¬ 
tration  is  approximately  zero  for  x2  +  y2  ^  r2.  then  a  circularly  symmetric  window  may  be 
appropriate.  Ii  would  work  very  well  if  the  Fourier  transform  data  had  a  circular  region  of 
support.  However,  spotlight  mode  SAR,  in  general,  does  not. 

4.1.2  Optimal  window*  for  SAR 

Windows  can  also  be  designed  which  are  optimal  in  a  given  sense.  In  [63].  Papoulis  listed 
several  optimization  criteria  such  as  maximum  energy  concentration,  specified  zero  crossings, 
minimum  energy  moment,  and  the  minimum  amplitude  moment.  The  maximum  energy  con¬ 
centration  criterion  generates  a  window  which  has  maximal  energy  in  a  given  range  ( — cr  to  <r) 
and  is  solved  via  the  prolate  spheroidal  wave  functions  (PSWF)..  Solutions,  both  continuous 
and  discrete,  of  the  PSWF  are  discussed  by  Slepian  [64.65.66].  The  specified  zero  crossings 
criterion  is  a  generalized  window  of  the  form 

w(t)  =*  k(  1  +  2a  cos(at))  pT(t)  (4.1) 

where 

-2sin(at)  / . 

2ar  +  sin(2ar) 


where  the  zero  crossings  are  specified  by  a..  The  Hamming  window  is  a  special  case  where  the 


Figure  4.1b  A  Circularly  Symmetric  Hamming  Window.: 


Figure  4.2a  A  Separable  Hanning  Window. 


Figure  4.2b 


A  Circularly  Symmetric  Hanning 


Window. 
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discrete  window  transform  has  zero  crossing  on  the  DFT  sample  points.  This  property  makes 
the  Hamming  window  useful  for  the  interpolation  evaluation  problem  since,  a  properly  placed 
and  sampled  Fourier  domain  Hamming  window  will  result  in  an  image  sampled  in  the  nulls  of 
the  image  domain  sine  function.  In  this  way.  the  interpolation  error  can  be  calculated  indepen¬ 
dently  of  windowing  effects.  Several  examples  of  circularly  windowed  data  records  are  given 
in  the  experimental  results. 

The  previous  sections  demonstrated  the  design  of  two-dimensional  windows  via  a  rather 
heuristic  algorithm  (  i.e..  take  a  good  one-dimensional  window  and  transform  it  into  a  two- 
dimensional  window).  However,  neither  method  takes  advantage  of  the  shape  of  the  region  of 
known  data.  O'Brien  [67]  and  Staehling  [68]  examine  the  question  of  optimizing  the  design  of 
windows  for  irregularly  shaped  data  regions.  Since  the  spotlight  mode  SAR  data  lie  on  a 
toroidallv  shaped  region,  it  is  worthwhile  to  examine  how  these  windowing  strategies  can  be 
used  prior  to  the  FFT  operation. 

O'Brien  investigated  the  use  of  an  energy  maximization  criterion  to  generate  windows 
which  can  be  used  on  spatially  limited  images  of  irregular  shapes,  e.g..  an  L  shape,  a  duel  cone, 
and  a  circle.  The  solution  is  based  on  an  iterative  application  of  the  time  and  frequency  limit¬ 
ing  operations  found  in  the  Papoulis  extrapolation  procedure  [69].  The  window  design  is  for¬ 
mulated  as  an  eigenvalue  problem. 

A  f(n.m)  »  B  T  f(n.m)  (4.3) 

where 

B  *  W-‘BW  is  the  bandlimiled  operator. 

T  is  the  spatial  truncation  operator. 

B  is  the  frequency-limited  operator,  and 

W  is  the  DFT  matrix. 


O'Brien  goes  on  to  demonstrate  the  convergence  of  the  algorithm,  and  gives  several  examples 
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which  show  the  types  of  windows  designed  for  the  simple  geometric  shapes  given  above.. 

Staehling  used  the  Hamming  window  applied  in  directions  orthogonal  to  the  edges  of  the 
given  data  set  (Fig.  4.3).  The  individual  window  arcs  were  expanded  to  meet  the  edge  of  the 
kno'*'r.  data  region.  His  results  demonstrated  better  sidelobe  control  than  the  separable  or  cir¬ 
cular  case.  This  may  be  expected  since  the  windowed  data  record  becomes  a  function  which 
gradually  approaches  zero. 

The  limitation  of  these  algorithms  is  due  to  the  lack  of  knowledge  of  the  shape  of  a  gen¬ 
eral  data  set.  Also,  the  algorithms  are  numerical  optimization  procedures.  Specifically,  the  data 
region  shape  must  be  known  before  design  or  application  of  the  adaptive  window  design.  If 
this  a  priori  information  is  unavailable,  then  the  separable  window  (somehow  optimized  for  the 
full  region)  is  the  better  and  (computationally)  simpler  alternative.  In  the  case  of  the  toroidal 


Figure  4.3  Application  of  the  Hamming  window  to  Orthogonal  Edges. 


data  set.  the  previous  strategies  could  be  used  only  if  the  toroid  were  sampled  on  a  rectangular 
raster  prior  to  windowing  with  zeros  surrounding  the  data  record  (Fig.  4.4).  This,  however,  is 
precisely  the  problem  faced  in  the  interpolator.  After  the  interpolation  step,  the  window  boun¬ 
dary  is  a  moot  point  since  the  data  are  now  on  a  rectangular  raster  and  standard  windows  can 
be  applied.  Because  the  Hamming  window  allowed  the  placement  of  zeros  at  the  image  sine 
zeros,  it  was  used  after  the  interpolation  stage  for  the  sidelobe  reduction. 


Figure  4.4  The  Toroidal  SAR  Data  Set  on  a  Rectangular  Raster, 


4.2  The  Limited  Data  Record  and  Extrapolation 

This  section  briefly  discusses  the  effects  that  a  limited  data  set  in  the  Fourier  domain  has 
on  the  image  reconstruction.  This  is  a  ubiquitous  problem  for  signal  processing  and  is  covered  a 
great  deal  in  the  area  of  extrapolation..  The  data  set  that  we  are  given  extends  over  only  a  very 
small  part  of  the  Fourier  domain  -  a  piece  of  a  torus.  Outside  this  region,  the  data  are  non-zero 
(for  spatially  limited  images,  it  cannot  be  zero),  but  in  interpolating  from  one  grid  to  another, 
it  is  often  assumed  to  be  so.  Bandlimited  extrapolation  of  the  known  signal  set  would  seem  to 
be  worthwhile  in  SAR  for  two  reasons:  1 )  the  resolution  of  the  system  is  proportional  to  the 
size  of  the  Fourier  piece  so  that  by  extending  the  data  set.  we  can  increase  the  system  resolu¬ 
tion:  2)  when  interpolating  from  one  grid  to  another,  the  high-order  interpolators  will  not  "fall 
off"  the  sample  grid  near  the  edges.  In  practice.  (2)  is  not  as  much  a  problem  as  would  appear, 
because  the  size  of  a  real  system  raster  is  several  orders  of  magnitude  larger  than  the  interpola¬ 
tor  kernel  ana  errors  would  only  occur  near  a  few  corners  of  the  grid  (Fig.  4.5).  O'Brien  [67] 
briefly  examined  the  possibility  of  extrapolating  the  Fourier  data  set  beyond  the  recorded 
region,  but  concluded  with  results  consistent  with  Abend  [70],  who  notes  that  bandlimited 
extrapolation  has  only  limited  use  in  spectral  estimation. 

4.3  Spurious  Targets 

It  is  well-known  that  interpolation  errors  in  the  Fourier  domain  lead  to  a  noisy  response 
in  the  spatial  domain.  The  very  nature  of  the  transform  operation  predicts  a  global  error  in  the 
spatial  domain  for  singular  errors  in  the  Fourier  domain.  The  following  sections  review  and 
attempt  to  qualitatively  characterize  this  interpolation  error  effect  for  tomography  and  SAR.. 

4.3.1  Munson-O'Brien  analysis  In  [67]  O’Brien  examined  the  effect  of  interpolation  in  one 
domain  on  the  transform  domain,  specifically  looking  at  sample  rate  changes  during  the  ;nler- 
polation  stage.  The  results  help  explain  the  presence  of  spurious  targets  in  the  output  image. 


Figure  4.5  Areas  where  the  Interpolation  Kernel  Falls  Off  the  Edge. 


In[l6],  Munson  clarifies  O'Brien's  analysis  of  spurious  targets.  The  result  of[l6].  are  presented 
briefly  here  for  completeness. 

The  analysis  can  be  simplified  if  the  Cartesian-to-Cartesian  (C-C)  interpolation  problem  is 
reduced  to  one-dimension..  Though  the  SAR  interpolation  problem  is  polar-to-Cartesian  (P-C) 
rather  than  the  simpler  C-C  geometry,  the  analysis  can  be  used  to  predict  where  spurious  tar¬ 
gets  may  appear,  since  the  section  of  the  polar  grid  is  approximately  rectangular.  This  type  of 
approximation  may  not  be  useful  in  the  tomographic  setting  where  the  grids  have  markedly 
different  shapes. 


The  sample  rate  conversion  problem  can  be  viewed  as  an  analog  interpolator  followed  by 
a  sampler  with  the  new  period  T'.;  In  the  case  where  the  interpolation  is  performed  in  the 
Fourier  domain,  the  process  can  be  viewed  as  shown  in  Fig.  4.6,  with  the  input  sample  set 


85 


available  after  the  initial  transform  operation.  The  case  and  hat  notation  here  is  like  that 
defined  in  the  Chapter  2  discussion  on  the  projection  slice  theorem.  The  imaging  process  per¬ 
forms  the  Fourier  transform  of  f(t).  into  F(w).  F(o>)  is  then  sampled  at  times  nT  to  form 
F'(nT).  The  sample  rate  conversion  step  can  be  thought  of  as  a  two  step  process  of  analog  inter¬ 
polation  followed  by  a  resampling  operation.  Because  the  interpolation  step  is  inexact,  it  pro- 

A 

duces  the  continuous  function  F(o>)  which  contains  aliased  copies  of  f(t).  This  is  resampled  as 

A  A 

F(mT')  and  converted  back  to  the  time  (spatial)  domain  as  f(k)  via  the  FFT  operation.  We 
would  like  to  determine  the  relationship  between  f(k)  and  f(x).  Note  that  f(k)  is  a  sampled 
function,  and  the  original  function  f(x)  is  continuous.  The  continuous  interpolation  kernel  is 
h(t)  with  transform  H(o>).  It  will  lie  convenient  to  define  the  infinite  sampling  sequence 

on 

ST(x)  =  £  8(.\— kT)  which  is  used  in  both  the  time  and  frequency  domain.  It  is  important  to 

-A  1* 

remember  that  the  operations  normally  performed  in  the  time  domain  (sampling,  convolution, 
etc.)  are  now  being  done  in  the  frequency  domain,  but  the  duality  properly  of  Fourier  theory 
can  be  used  to  understand  what  is  happening  in  the  image  (time)  domain. 

The  output  of  the  initial  sampler  (after  the  Fourier  inversion  step)  is 

F(nT)  =  F(w)  •  Sx(<a).  This  has  the  effect  of  convolution  with  the  infinite  impulse  sequence 

S2w(t)  in  the  time  domain.  Aliasing  can  result  if  T  is  insufficiently  small.  The  sampled  Fourier 
T 

sequence  is  convolved  with  the  continuous  interpolation  kernel.  H(  omega  )  to  reconstruct  an 
approximation  to  F 


F(w)  =  £  F(nT)  H(g>  -  nT) 

n 

=  |F(o>)-ST(oi)}  *  H(oi) 


(4.4a) 

(4.4b) 


In  the  time  domain,  this  becomes 

{f(t)*S2ff(t)}  •  h(t) 

f 


(4.5) 


The  final  sampling  step  multiplies  the  reconstructed  function  with  the  infinite  sampling 
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W-1701 


Figure  4.6  Sample  Rate  Conversion  Model  for  Fourier  Data  Interpolation. 

function,  but  with  period  T 

{ F(w)Vw)}  *  H(«)  |  •  Sr(«)  (4.6) 

and  in  the  time  domain 

{f(t)*S^(t)}-h(t) 

T 

It  is  precisely  this  second  convolution  after  the  h(t)  windowing  which  accounts  for  spurious 
(aliased)  targets  to  appear  in  the  reconstructed  image.  It  was  shown  in  the  interpolation 

2tt 

chapter  that  the  interpolation  kernel  must  have  a  transform.  h(t),  which  removes  the  -y- 

spaced  copies  of  the  spectrum  (In  this  case,  the  copies  are  of  the  image.)  If  h(U  were  bandUm- 
ited.  the  spurious  targets  would  not  appear.  This  would,  of  course,  create  an  infinite  length 
interpolation  kernel. 
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4.3.2  Analysis  for  the  nearest  neighbor  interpolator 

Stark  analyzed  the  nearest  neighbor  algorithm  in  terms  of  random  jitter  noise  in  the 
polar-rectangular  interpolation  stage  to  explain  streak  artifacts  and  high-frequency  noise 
induced  into  a  direcl-Fourier  reconstructed  tomographic  image  [22].  He  generates  a  noise 
model  of  the  positional  error  in  the  transform  domain  as  a  uniformly  distributed  variable 
between  -1/2  and  +1/2.  This  corresponds  to  the  position  error  induced  from  using  the  nearest 
neighbor  which  must  be  within  1/2  a  sample  point  away  (sample  period  T  normalized  to  1). 
This  positional  error  “noise"  is  uncorrelated  with  the  signal..  Stark  showed  a  relationship 
between  the  positional  error  noise,  the  Fourier  transform  of  the  original  signal  and  the  Fourier 
transform  of  the  interpolated  signal  which  demonstrates  that  the  positional  error  generates 
noise  which  is  correlated  with  the  signal.  He  conciudeed  that  the  correlated  noise  translates  into 
streaking  artifacts  in  the  spatial  domain,  as  well  as  producing  high-frequency  amplification. 
The  streak  artifacts  are  also  found  in  the  SAR  interpolation  problem  as  smearing  across  the 
reconstruction.. 

Some  insight  into  what  is  occurring  can  be  gained  by  looking  at  the  u.v  error  terms  of  each 
interpolated  point  for  the  nearest  neighbor  algorithm.  Let  E„  be  the  amount  of  error  in  the  u- 
direction  and  Ev  be  the  error  in  the  v-direction  (Fig.  4.7)  in  the  Fourier  plane.  In  general,  Eu 
and  Ev  will  both  be  functions  of  u  and  v.  Let  there  be  a  point  target  at  (u„  .  v„).  The  Fourier 
transform  is 


F 


8(u— u,„v— v„) 


F(U.V)  =  e-j(l'u°  +  Vv®:i 


(4.8) 


A  A  A  A 

Let  (U.V)  be  the  nearest  neighbor  to  the  polar  sample  (U.V).  Thus.  F(U.V)  is  used  instead  of 
the  true  F(U,.V\  Define  the  transform  relationship 

f(u.v)  =  f’-,{F(U,V)}  <=>  F(U.V)  =  F(U.V)  (4.9) 


Now.  expanding  F(U.V)  in  a  2-D  Taylor  series 
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»ut 


F(O.V)  =  F(U.V)  +  (U  -  L) 

eu 

+  SM*1  (V-V) 

dV 

+  higher  order  terms. 


(4.10) 


liu(U.V)  «  (0  -  U)  =  /•’2_D(eu(u.v)}  (4.11a) 

Ev(U.V)  =  ( V  -  V)  =  F2_D{ev(u.v)|  t4.1  lb) 

is  the  transform  of  the  positional  errors  for  each  coordinate  which  depends  on  (u,v).  The  par¬ 
tial  derivatives  of  F  are 


Thus.  (4.10)  becomes 

F(U.V)  *  = 


ai;(i;.v) 

df 

0FOJ.V) 

dV 


ju„  F(U.V) 
jv„  F(L'.V) 


F(U.V)  +  ju„  F(t'.V)  EU(U.V) 

+  jv„F(U.V)Ev(U.V) 


(4.12a) 

(4.12b) 


(4.13) 


where  the  higher  order  terms  have  been  dropped.  The  spatial  domain  equivalent  is 

f(u.v)  «  f(u.v)  +  juo  f(u.v)  **  eu(u,v)  (4.14) 

+  jv0  f(u.v)  **  ev(u.v) 

Equation  (4.14)  shows  that  the  point  target  position  and  2-D  error  surface  are  related  to  the 
image  reconstruction  though  a  2-D  convolu'.jon,  and  that  the  error  surface  is  proportional  to  the 
position  of  the  point  target  u„  and  v„. 


Now  as  a  special  case,  suppose  E„(U.V)  and  Ev(L.V)  are  2-D  periodic  with  fundamental  fre¬ 
quencies  (au.  0U)  and  (0U.  0V).  Then:; 
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EU(U,V)  =  £1  a*  ejtlcf"u  +  k“vV> 

i  k 

Ev(U.V)  =  Z  I  bIk  ew’  +  k^v) 

i  k 

and 


(4.15a) 

(4.15b) 


f(u.v)  3:5  f(u.v)  +  juii  II  aik  f(u  +  iau.v  +  kav)  (4  16) 

i  k 

+  jv0  II  bik  f(u+i0u,v+kafv) 

i  k 

Equation  (4.16)  represents  the  special  case  when  the  error  surfaces  in  Fourier  space  are 
periodic.  If  the  sampling  frequencies  in  the  input  sampling  array  are  approximately  constant, 
as  in  the  case  where  only  sub-patches  are  used  in  the  reconstruction,  then  the  interpolation  will 
be  from  a  near-rectangular  grid  to  an  exact  rectangular  grid,  which  is  perhaps  rotated.  The 
error  surfaces  in  this  case  are  periodic  triangle  waveforms  with  an  orientation  that  determines 
a  and  0  (Figs.  4.7a  and  4.7b).  By  knowing  the  relative  geometric  orientations  of  the  grids,  the 
resultant  smear  pattern  may  be  predicted. 

Equation  (4.16)  also  shows  that  the  smear  magnitude  is  proportional  to  the  displacement 
of  the  point  target  from  the  origin.  Points  farther  away  from  the  patch  center  will  cause 
greater  smear.  This  is  similai  to  the  effect  caused  by  the  quadratic  phase  term  in  the  SAR  phase 
equations  (Chapter  2). 

The  error  surface  directional  vectors  defined  as  a  *  (au.av)  and  $  *  (0U.0V)  describe  the 
smearing  direction  and  depend  on  the  error  surface  orientation.  The  input  polar  grid  displays  a 
spatially  varying  period  and  orientation  which  is  slowly  modulated  with  respect  to  the  rec¬ 
tangular  grid.  Although  these  error  surfaces  may  not  exhibit  2-D  periodic  behavior  in  the  gen¬ 
eral  case,  for  the  near  rectangular  SAR  geometry,  and  rectangular  to  rotated-rectangular  grid 
geometries,  they  are  approximately  periodic. 
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4.4  Alternative  Sampling/PRF  Strategies 

The  previous  sections  have  dealt  primarily  with  the  standard  polar  input  raster. 
Modifications  to  the  sampling  and/or  radar  pulsing  mechanism  allows  other  sampling  rasters  to 
be  used  with  greater  efficiency  and  reduction  of  error.  This  concept  was  presented  in  a  tomo¬ 
graphic  setting  by  Mersereau  and  Oppenheim  [71]  through  the  use  of  concentric  square  sam¬ 
pling.  The  resulting  image  reconstruction  is  significantly  improved  because  interpolation  is 
only  done  in  one  direction.  Simplifying  the  interpolation  geometry  to  one-dimension  also 
reduces  the  computation  of  the  interpolation  stage,  both  in  locating  the  grid  points  and  per¬ 
forming  the  filtering/sample  rate  changes.  Of  course,  this  type  of  sampling  strategy  requires  a 
much  more  sophisticated  sampler.  In  SAR.  alternate  grids  may  be  designed  to  take  advantage 
of  a  smart  sampling  device.as  well  as  a  programmable  pulse  repetition  frequency  (PRF)  genera¬ 
tor.  The  keystone  grid  is  rather  well-known,  but  more  as  an  intermediate  interpolation  grid 
rather  than  an  input  raster.  Another,  is  a  polar-like  grid  with  equi-PRF  spacing.  The  range 
samples  are  uniformly  sampled  for  each  return  signal,  and  the  pulses  are  transmitted  at  a  uni¬ 
form  rate  (the  PRF  is  constant)  as  the  radar  moves  past  the  terrain  at  a  constant  velocity.  This 
results  in  a  range  lines  which  are  unequally  space  in  9.,  However,  the  intermediate  (keystone) 
grid  produced  by  the  separable  algorithms  will  have  equally  spaced  samples  in  the  azimuthal 
direction. 

4.5  Keystone  with  Smart  PRF 

If  the  incoming  signal  is  sampled  on  the  keystone  grid  shown  in  Fig.  3.6,  then  the  complex 
two-dimensional  interpolation  reduces  to  one-dimension.  This  greatly  reduces  interpolation 
error  and  compulation  lime  as  shown  by  the  computer  simulations  of  Chapter  5.  It  also  allows 
for  a  novel  interpolation-Fourier  transform  operation  in  one  step  -  the  Chirp  7-transform. 


4.5.1  The  chirp  z-transform  algorithm 

An  N-poinl  DFT  of  a  sequence  can  be  thought  of  as  sampling  the  z-transform  of  a 
sequence  uniformly  around  the  unit-circle  in  the  z-plane.  The  chirp  z-lransform  is  widely 
known  as  a  generalized  DFT  algorithm  that  computes  M  samples  of  the  z-transform  along  an 
arbitrary  spiral  contour  in  the  z-plane.  Once  more.  N  and  M  need  not  be  the  same  as  in  the 
DFT.  nor  do  they  have  to  be  highly  composite  for  implementation  with  an  FFT.  The  chirp  z- 
transform  is  based  on  using  the  FFT  to  perform  fast  convolutions  of  the  input  sequence  with  a 
chirp,  or  frequency  ramped,  signal.  The  algorithm  is  explained  in  its  entirety  by  Rabiner  et  al. 
in  [72]  of  which  a  distilled  version  is  found  in  [73].  It  is  used  in  pole  enhancement,  narrow- 
band  frequency  analysis,  bandlimited  interpolation,  and  arbitrary  radix  DFT  computations. 
Although  the  most  general  form  of  the  chirp  z-transform  allows  for  computation  on  an  arbi¬ 
trary  spiral  contour  in  the  z  plane,  this  application  uses  it  to  simply  perform  an  M  point  D1T 
with  an  N  point  sequence.  M  and  N  being  different.  This  means  that  the  contour  reduces  to  the 
unit  circle  in  the  z-plane  with  the  output  sequence  beginning  at  a  non-zero  phase  angle. 

The  M  point  chirp  z-transform  of  the  N  point  sequence  x(n)  is  given  by 

X(m)  »  Nf  A”nWnm  (4.17) 

n®o 

where  A  *  Aoe*2*^'  specifies  the  output  start  angle  and  radius,  and  W  *  W(tei2w*°  specifies  the 
contour  curvature,  and  output  sample  spacings  as  shown  in  Fig.  4.8.  In  the  SAR  case,  the  out¬ 
put  contour  is  still  on  the  unit  circle,  although  it  may  be  located  at  the  arbitrary  angle  2irB„ 
and  extend  to  2ir(M— l)0ft  where  is  the  output  sampling  angle  (  ;ig.  4.9). 

Through  some  algebraic  manipulation,  (4.17)  can  be  reduced  to  a  3-step  process  of 

(1)  Ramping  (multiplying  by  a  frequency  swept  sequence)  the  !\  point  input  sequence  by 
A~"  W“J/2. 

(2)  Convolving  this  result  with  W-0^2., 
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(3)  Deramping  the  convolved  sequence  with  VV,,2/2. 

This  is  illustrated  in  Fig.  4.10. 

The  normal  convolution  can  be  done  with  an  FFT  if  the  sequences  are  properly  padded 
with  zeros  to  avoid  wrap  around  (circular  convolution).  Unfortunately,  the  zero  padding 
increases  the  FFT  size  to  N'  which  is  the  first  power  of  «.wo  greater  than  N  +  M.  This  can  lead 
to  much  higher  computation  time  for  small  N.  and  make  the  chirp-z  algorithm  complexity 
disproportionately  high  compared  to  the  other  interpolation  algorithms.  It  must  be  remem¬ 
bered.  however,  that  the  chirp-z  approach  folds  the  Fourier  domain  inversion  step  into  the 
interpolation  stage  so  that  comparisons  must  be  made  with  the  2-D 1FFT  step  added  to  the  com¬ 
plexity  measure.  This  becomes  less  of  a  problem  as  N  is  increased,  since  the  FFT  time  is  order 
NlogjX. 


h(n)-\r*n,/2 


rigure  4.10  The  Chirp  7-transform  as  a  Pre-  and  Post-Multiplied  Convolution. 
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In  the  SAR  geometry,  the  Keystone  sampling  and  smart  PRF  place  data  points  on  straight 
azimuth  lines  with  sample  points  equi-spaced  on  any  particular  azimuth  line  (though  the  spac¬ 
ing  changes  from  azimuth  line  to  azimuth  line).  The  chirp-z  algorithm  is  applied  to  each  key¬ 
stone  azimuth  line  placing  the  output  points  on  the  rectangular  grid.  A  standard  FFT  is  per¬ 
formed  in  the  orthogonal  direction.  Here.  A  =  exp(j0„)  and  W  =  exp(— #o).  where  0O  is  the 
angle  of  the  first  spectral  sample  and  0„  defines  the  output  sample  spacing. 

The  number  of  input  points  N  for  the  chirp-z  algorithm  is  arbitrary,  and  the  number  of 
output  points  \l  is  equal  to  the  rectangular  raster  size  -  1024  for  the  total  image.  64  for  the 
computer  simulations.  Fach  azimuthal  input  line  is  windowed  (Hamming)  prior  to  the  chirp 
z-transform.  Note  that  as  the  azimuth  line  index  increases  with  range,  the  input  sample  spacing 
increases  and  the  window  size  will  increase.  This  is  a  modification  of  the  warped  windows  dis¬ 
cussed  by  Slaeling  and  O’Brien.  The  range  line  window  is  applied  prior  to  the  secondary  FIT 
stage. 

The  input  data  set  lies  on  a  trapezoidal  raster  and  the  output  of  each  chirp  z-transform 
must  be  correctly  phased  between  azimuth  lines.  This  can  be  done  by  multiplying  each  output 
point  by  a  linear  phase  term  that  changes  from  azimuth  <  ne  to  azimuth  line 

expljv  (u  tanOI 

Because  the  chirp  z-transform  is  a  forward  transform  (which  is  followed  by  a  forward 
FFT)  and  the  data  set  is  already  in  the  Fourier  domain,  the  resulting  image  must  be  time 
reversed,  i.e..  rotated  180  *.  to  correct  the  effect  of  a  double  forward  transform  since 


FT 


FTlf(x.y)} 


f(-x.-y) 


4.6  Polar  Format  with  Equi-PRF 

A  slight  modification  to  the  PRF  will  place  the  samples  on  an  equi-PRF  grid  which  appears 
like  a  polar  grid,  but  with  angular  increments  that  are  not  constant  (Fig.  4.11).  This  has  the 
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advantage  that  the  separable  algorithms  in  Chapter  3  which  generate  an  intermediate  keystone 
grid  will  have  equi-spaced  sample  points  along  the  azimuth  lines  (just  as  the  smart-PRF  key¬ 
stone  grid  in  the  previous  section).  This  simplifies  the  interpolation  step  in  the  second  stage  of 
the  separable  algorithms  (weighted  sine,  spline)  and  thus  reduces  computation. 


Figure  4.11  Equi-PRF  Polar  Grid. 
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CHAPTER  5 

EXPERIMENTAL  EVALUATIONS 

The  previous  chapters  presented  the  theory  of  interpolation  for  signal  reconstruction  in 
both  spatial  and  frequency  domains.  This  chapter  applies  this  developed  theory  to  the  original 
problem  of  image  reconstruction  in  synthetic  aperture  radar. 

The  following  sections  describe  the  approach  taken  for  the  computer  experiments,  the 
simplification  to  the  SAR  model  used,  and  the  results  of  the  various  interpolators.  A  novel 
measure  of  image  reconstruction  is  also  presented  here.  Because  the  number  of  parameters  is 
large  (Fourier  section  of  interest,  point  ta.get  location,  data  record  window,  interpolator  type, 
and  interpolator  order)  only  a  subset  of  image  plots  :>  ,iven.  A  more  complete  set  of  interpola 
tor  performance  results  are  given  in  tabular  form  (Tables  C.l  to  C.12)  in  Appendix  C. 

5.1  Interpolation  Model 

In  the  evaluations,  we  assume  a  look  angle.  0„,M  of  3*.  This  is  the  angle  over  which  all  of 
the  data  are  collected.  The  rectangular  array  to  be  interpolated  onto  is  inscribed  within  the  two 
arcs  Rmm  and  Rm,x.  which  are  proportional  to  the  chirp  bandwidth  defined  in  Chapter  2.  The 
input  polar  sample  array  is  assumed  to  have  a  size  of  1024  by  1024  points  with  1024  equally 
spaced  samples  along  each  range  line  R„,,„  to  Rnwx  and  1024  samples  along  arcs  from  — dmax/2  to 

A  program  to  generate  SAR  data  from  the  equations  of  Chapter  2  was  written,  allowing 
each  of  the  many  geometric  and  system  parameters  to  be  adjusted.  While  (2.26).  represented 
an  exact  point  target  response,  the  smearing  caused  by  the  quadratic  phase  term  obscured  the 
effects  of  the  interpolator,  and  thus  were  abandoned  for  a  newer,  simpler  model:  the  Fourier 
transform  of  a  spatially  offset  impulse  function.  A  single,  'deal  point  target  was  placed  at 
(x0.Vo)  and  the  corresponding  Fourier  transform.  e_jU  +  Vy«>i  was  sampled  on  the  uniform 
polar  grid.  Allowable  target  positions  which  did  not  produce  aliasing  ranged  from  -32  to  +31 
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in  the  spatial  domain.  We  used  (-23,24)  as  the  point  target  position  in  the  evaluations 
presented  here.  This  stressed  the  interpolator  by  including  high-frequency  components  while 
allowing  the  filter  response  to  have  a  finite-width  transition  band.  High  frequencies  stress  the 
interpolator  because  a  higher  degree  of  variation  is  more  difficult  to  reproduce  with  simple  ker¬ 
nels. 

Because  our  computer  resources  did  not  permit  processing  such  a  large  array,  we  chose  to 
examine  interpolator  effects  on  various  64  x  64  subarrays  of  the  full  1024  by  1024  grid.  Each 
64  x  64  subarray  is  numbered  from  1  to  16  along  t he  U-axis.  and  from  1  to  8  along  the  V-axis. 
For  example,  subarrays  (2.8)  and  (16,1)  are  shown  in  Fig.  5.1.  Symmetry  about  the  U-axis 
eliminates  the  need  to  study  the  interpolator  for  v  <  0.  The  subarrays  (2.8),  (10.5),  (16.1)  and 
(16.8)  were  studied,  since  they  represented  various  sections  of  the  full  array  that  would  have 
unique  properties,  i.e..  rotation,  maximal  sample  rate  change,  average  performance,  etc.  The 
input  data  region  cf  each  subarray  was  extended  in  both  range  and  azimuth  by  15  sample 
points  in  each  direction  (making  the  available  data  94  by  94  samples)  so  that  the  higher  order 
interpolators  did  not  "fall  off"  the  edge  of  the  input.  .Since  the  original  1024  by  1024  grid  had 
only  a  very  small  set  of  points  where  this  occurs,  this  data  extension  was  justifiable  for  the 
subarrays. 

5.2  A  Novel  Fijure-of -Merit 

Evaluation  of  the  interpolators  is  a  difficult  task,  especially  since  radar  data  interpretation 
tends  to  be  subjective.  The  figure  of  merit  that  we  have  used  is  called  the  Multiplicative  Noise 
Ratio  (MNR).  which  is  very  similar  to  the  more  commonly  known  "integrated  sidelobe  ratio” 
frequently  used  for  one-dimensional  radar  evaluation  [l].  Because  the  rectangular  grid  is 
weighted  with  a  separable  Hamming  window,  an  ideal  point  target  is  spread  3  samples  in  each 
dimension  U  and  V.:  In  the  following  experiments,  the  mainlobe  of  the  response  is  defined  as 
the  5x5  square  centered  on  the  peak  response.  The  MNR  for  the  reconstruction  is  then  defined 
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N  x  N  subarray  (16,1) 

Figure  5.1  Geometry  of  Fourier  Domain  Subarray  Locations. 
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Figure  5.2  Exact  Reconstruction  for  a  Target  at  (-23.24), 


MNR  =  10-log 


^(magnitude  of  the  points  outside  mainlobe)2 
^(magnitude  of  the  points  inside  mainlobe)2 

Figure  5.2  shows  a  perfectly  reconstructed  ideal  point  target  at  (-23.0.24.0)  with  a  calculated 
MNR  of  -48.13  dB.  It  was  produced  by  generating  target  data  directly  on  a  rectangular  grid 
and  bypassing  the  interpolating  stage  during  reconstruction.  The  floor  for  all  point  target 
response  plots  has  been  set  to  -60  dB  relative  to  the  peak  value.  Note,  in  Fig.  5.2.  the  effect  of 
the  Hamming  window  on  the  ideal  response.  The  point  target  has  widened  and  the  rest  of  the 
samples  fall  on  the  nulls  of  the  sine  (the  response  for  a  point  target  on  a  limited  record).  If  the 
target  is  moved  to  a  position  halfway  between  sample  coordinates,  the  sine  response  is  then 
sampled  on  the  peaks  of  the  sidelobes.  resulting  in  a  badly  distorted  response  and  an  MNR  of 
-28.5  (Fig.  5.3).  This  demonstrates  that  the  MNR  calculations  are  valid  only  for  targets  which 
are  positioned  directly  on  the  sample  points.  If  the  interpolator  caused  the  target  to  move  only 
slightly,  the  sidelobes  become  visible  and  the  MNR  increases  misleadingly.  It  is  interesting  to 
note  that  the  MNR  can  produce  values  which  apparently  differ  from  a  subjective  view.  Some¬ 
times.  one  reconstruction  will  look  better  than  another,  but  the  MNR  value  is  worse  (more  posi¬ 
tive).  This  is  usually  the  result  of  a  very  wide  target  response  in  which  the  samples  directly 
outside  the  true  5  by  5  mainlobe  are  summed  in  the  denominator,  causing  an  increase  (more 
positive)  in  the  MNR.  The  results  are  thus  obtained  as  follows: 

(1)  Generate  a  sampled  polar  array  by  sampling  the  Fourier  Transform  of  the  ideal  target. 

(2)  Interpolate  to  the  rectangular  grid  of  interest  (1.1)  to  (16.8). 

(3)  Window  the  resulting  64  X  64  rectangular  data  subarray  with  a  separable  Hamming  win¬ 
dow. 

(4)  Calculate  the  IFFT. 


i()2 


(5)  Calculate  the  MNR  of  the  spatial  domain  image. 

(6)  Display  the  log10  of  the  squared  magnitude  of  the  reconstruction  (the  power)  with  a  floor 

of  -60  dB. 

The  program  evaluate  calculates  the  image  MNR  and  locates  the  largest  (magnitude) 
twenty  points  on  the  grid  to  determine  the  positions  of  the  spurious  peaks  and  false  targets. 
The  dB  levels  of  these  points  relative  to  the  detected  peak  are  also  calculated. 

5.3  Algorithm  Complexity 

This  work  focused  on  improving  the  SAR  reconstruction  quality,  while  reducing  the 
interpolation  processing  time.  i.e..  reducing  the  computational  resources  required  to  perform  the 
algorithm.  The  required  resources  are  specified  by  what  we  have  termed  "complexity"  -  a 
measure  of  how  complex  (in  terms  of  mathematical  operatations)  the  algorithm  is  to  perform. 
Algorithm  complexity  is  a  difficult  issue,  because  the  algorithms  presented  do  not  offer  an  order 
of  magnitude  increase  in  processing  speed,  but  typically  improve  speed  by  a  factor  in  the  range 
of  2  to  10.  Also,  the  results  are  dependent  on  the  order  of  the  interpolator. 

Crochiere  and  Rabiner{35]  use  the  notion  of  Multiply-ADditionS  or  MADS/sec  when 
optimizing  the  decimation/interpolation  stages  for  a  sample  rate  converter.  This  was  used 
because  the  final  design  was  a  simple  time-invariant  FIR  filter  with  real  coefficients  which  has 
only  multiplies  and  adds.  Neither  complex  arithmetic  or  divisions  are  required,  nor  transcen¬ 
dental  function  evaluation  (sine,  square  root.  etc.).  They  also  make  the  assumption  that  an  add 
takes  the  same  amount  of  time  as  a  multiplication.  This  may  not  be  totally  unfounded,  as  in 
the  case  when  all  arithmetic  (necessarily  integer)  is  done  with  a  look-up  table  and  requires  only 
one  clock  cycle.  Real  multiplies,  however,  typically  require  more  time  than  real  adds,  and  real 
divide  operation  is  more  expensive  than  a  real  multiply  (although  no  divides  are  performed  in 
the  simple  ID  FIR  filters). 
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Here,  the  algorithm  complexity  is  calculated  by  computing  the  number  of  operations  of  a 
given  type  that  occur  in  the  interpolation  software.  Though  this  measure  would  seem  to  be 
based  on  the  particular  programming  practice,  it  takes  into  account  many  of  the  operations  nor¬ 
mally  ignored  by  "order  of  complexity  analyses,"  such  as  calculation  of  the  coordinate  spatial 
position  and  the  evaluation  of  trigonometric  functions.  The  complexity  was  determined  by 
diagramming  program  loops  and  generic  operations  and  identifying  function  calls.  Tri¬ 
gonometric  and  transcendental  function  calls  were  given  their  own  unit  of  cost  measure  so  that 
the  complexity  would  not  be  skewed  by  a  particular  implementation,  such  as  look-up  table, 
power  series  expansion,  or  a  table-interpolation  scheme.  Thus,  the  cost  to  evaluate  the  sin 
function  is  represented  as  C,rig.  the  cost  to  compute  a  square  root  is  CS(|rt.  the  cost  to  compute  a 
complex  exponential  is  Cw,„  and  the  cost  to  generate  a  sine  value  is  Csmc.  The  cost  of  a  real  add 
or  subtract  is  C,/,  and  the  cost  of  a  real  multiply  or  divide  is  Cm/<1.  Although  it  may  seem  that  • 
this  is  a  crude  estimate  of  algorithm  complexity,  it  has  been  found  through  our  experience  to  be 
reasonable.  The  complexity  calculation  was  based  on  the  inlet polaior  alone,  and  did  not 
include  the  windowing,  file  reading/writing,  or  post  processing  limes  (array  permuting,  magni¬ 
tude  detection,  etc.).  The  FFT  stage  has  been  included  in  the  complexity  analyses  because  the 
chirp-2  algorithm  performs  both  interpolation  and  FFT  as  one.  inseparable  operation. 

The  computation  time  for  the  interpolation  set  was  measured  for  each  algorithm  with 
varying  orders  and  different  input  grids.  Because  the  VAX  CPU  timing  facility  is  system  load 
dependent  (the  measured  epu  time  is  a  function  of  the  load  of  the  system),  an  average  run  time 
value  was  computed  over  all  of  the  runs  for  a  particular  grid  and  interpolator.  This  value  is 
used  to  plot  the  algorithm  times  and  to  compute  another  performance  measure  -  the 
1MNR1/CPU  ratio.  This  ratio  is  useful  to  describe  the  quality/cost  of  an  interpolator.  Higher 
ratios  indicate  a  better  reconstruction  for  the  amount  of  computation  used.  It  is  only  useful  to 
compare  interpolators  within  a  given  subarray,  since  the  MNR  ratios  vary  dramatically 
between  subarrays  with  similar  order  interpolators. 
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5.4  Polar  Grid 

The  polar  grid  input  array  represented  the  most  well  known  data  geometry.  The  Fourier 
piece  (1,1)  is  the  subarray  which  is  closest  to  the  output  grid  in  both  sample  spacing  and  orien¬ 
tation.  This  leads  one  to  believe  that  the  nearest  neighbor  interpolator  should  perform  quite 
well  here.  In  fact,  it  does  very  well,  as  shown  in  Fig.  5.4  with  an  MNR  of  -42.12.  This  is  only  6 
dB  away  from  the  ideal  response. 

In  subarray  patch  (2.8).  the  nearest  neighbor  dots  not  perform  nearly  as  well  (Fig.  5.5). 
With  an  MNR  of  -5.64.  the  reconstruction  is  extremely  noisy  and  contains  many  streaky 
artifacts.  It  is  interesting  to  note  that  the  nearest  neighbor  produces  a  streak  pattern  which  is 
similar  to  the  artifacts  described  by  Stark  [22]  in  his  nearest  neighbor  analysis.  It  can  also  be 
explained  in  terms  of  the  nearest  neighbor  analysis  presented  in  Chapter  4.  The  rotated  grid  in 
(2.8)  generates  an  almost  periodic  positional  error  surface  which  has  an  orientation  related  to 


TABLE  5.1  Evaluation  Results  for  Fourier  Piece  (2.8)  on  Polar  Grid. 


Target 

Position 

Inlerp. 

Parameter 

MNR  (db) 

CPU  lime 
(seconds) 

1 MNR/CPU  1 

-23.24 

NN 

-5.64 

3.28 

1.72 

-23.24 

G-ID 

(0.1) 

-14.56 

19.82 

0.73 

-23.24 

G-ID 

(0.2) 

-18.41 

5.75 

3.20 

-23.24 

G-ID 

(0.3) 

-17.00 

22.57 

0.75 

-23.24 

wsinc 

2 

-18.92 

17.7 

1.07 

-23.24 

wsinc 

4 

-24.50 

28.4 

0.86 

-23.24 

wsinc 

6 

41.3 

0.74 

-23.24 

wsinc 

8 

-36.53 

53.0 

0.69 

-23.24 

wsinc 

10 

-42.65 

61.7 

0.69 

-23.24 

wsinc 

12 

-46.88 

70.6 

0.66 

-23.24 

wsinc 

14 

-48.06 

82.7 

0.47 

-23.24 

wsinc 

16 

-48.24 

94.7 

0.51 

-23.24 

B-spline 

Knrrrni 

mam 

■BSH 

-23.24 

B-spline 

IB 

I 

-23.24 

B-spline 

-0.75 

17.22 

1.60 

-23.24 

B-spline 

-1.00 

(EES 

17.22 

1.77 

-23.24 

■aswaransw 

-32.41 

27.7 

1.17 
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POLAR:  NEAREST  NEIGHBOR 


WINOOW*  HAMMING  1 

MNR  UN  OB)  :  -S. 64366 

FOURIER*  <  2.  81  TARGETt  -23.00.  24.00  MAGITUOE*  1.0 


Figure  5.5  Nearest  Neighbor  Interpolator  at  (2.8). 
POLAR:  NEAREST  NEIGHBOR 


MNR  UN  08!  *  -5.31121 
INTERPOLATOR:  NEAREST  NEIGH8CR 


FOURIER*  U0.  SI  TARGET;  -23.00.  24.00  HAGITUOt:  1.3 


Figure  5.6  Nearest  Neighbor  Interpolator  at  ( 10,5). 
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the  relative  input/output  grid  angles.  The  same  type  of  smearing  is  presen’  when  the  target  is 
moved  to  (1,-15)  (Fig.  5.7).  The  MNR  is  a  little  better  at  -10.63.  because  the  Fourier  data  is 
more  slowly  varying,  but  the  reconstruction  noise  obscures  the  image.  Finally,  in  the  (16.1) 
subarray,  the  nearest  neighbor  interpolator  generates  many  spurious  targets  due  to  the  change 
in  sample  rate  (Fig.  5.8).  The  azimuthal  rate  change  at  this  position  in  space  is  approximately 
1.05,  and  the  Munson-O'Brien  equations  [16]  predict  spurious  targets  along  the  axis  of  rate 
change  (the  predicted  position  is  27.2,  the  measured  position  is  27)  with  magnitudes  propor¬ 
tional  to  the  interpolator  frequency  response  beyond  the  cutoff  frequency.  Since  the  nearest 
neighbor  interpolator  transform  is  a  sine,  the  sidelobes  are  very  high,  and  thus,  produce  the 
great  numbers  of  high  spurious  targets.  It  is  difficult  to  measure  the  angular  orientation  of 
these  error  surfaces  (they  are  not  at  the  polar  grid  angle)  but  the  angled  smearing  is  quite  clear 
(Fig.  5.5).  The  same  type  of  smearing  is  present  in  the  nearest  neighbor  reconstruction  of 
subarray  (10.5):  however,  it  is  in  a  different  direction  (Fig.  5.6). 

Figure  5.7  also  shows  that  the  use  of  the  FFT  as  an  approximation  to  the  Fourier 
transform  has  produced  a  reconstruction  which  is  apparently  one  period  of  a  2D  periodic  signal. 
The  other  reconstructions  reflect  the  same  conclusion.  This  is.  of  course,  expected,  since  the 
finite  record  sampling  operation  implicitly  creates  a  periodic  signal. 

The  inverse  distance  (ID)  interpolator  performed  better  than  tht  nearest  neighbor  (NN)  in 
all  subarrays,  based  on  the  MNR  values.  In  (2,8)  (Fig.  5.9),  there  is  a  marked  decline  in  the 
MNR:  -14.56  dB  (ID)  compared  to  -5.64  dB  (NN).  Compare  this  with  the  inverse  distance 
squared  (ID2)  reconstruction  of  Fig.  5.10.  The  MNR  of  -18.42  suggests  that  the  ID2  is  slightly 
better  than  ID.  and  furthermore,  the  interpolation  time  is  3.45  times  faster.  This  is  a  significant 
improvement  overall  as  seen  by  the  IMNRI/CPU  ratio  in  Table  5.1. 

If  the  generalized  inverse  distance  interpolator  is  used  by  including  the  next  16  nearest 
data  points,  both  ID  and  ID2  reconstructions  are  poorer  (Figs.  5.11  and  5.12,  and  Tbl.  5.1). 
Increasing  the  order  of  the  inverse  distance  interpolator  beyond  2  did  not  seem  to  improve  the 
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POLAR:  NEAREST  NEIGHBOR 


HI NOON:  HAMMING 
MNR  (IN  06)  :  -10.62914 
NEAREST  NEIGHBOR 

FOURIER:  (  2.  8)  TARGET:  1.00.-15.00 


MAGITUOE*  1.0 


Figure  5.7  Nearest  Neighbor  Interpolator  at  (2.8).  target  at  (1  .-15). 


POLAR:  NEAREST  NEIGHBOR 


HINOOW:  HAMMING 
MNR  (IN  08)  :  -1.37431 
INTERPOLATOR:  NEAREST  NEIGHBOR 

COURIER:  (IS.  1)  TARGET:  -23.00.  24.30  MAGITUDE:  1.3 

Figure  5.8  Nearest  Neighbor  Interpolator  at  ( 16.1 ). 
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response  (see  Table  5.1 ).  This  is  verified  by  comparing  the  respective  MNRs. 

The  ID  and  ID2  interpolators  performed  only  3  dB  better  than  NN  in  the  (16,1)  subarray 
(Figs.  5.13  and  5.14).  Subjectively,  the  NN  reconstruction  at  (16.1)  (see  Fig.  5.8)  is  extremely 
poor  compared  to  ID  and  ID2,  and  yet  the  MNR  differ  by  only  a  few  dB.  Part  of  the  reason  for 
this  lies  in  the  fact  the  images  are  plotted  with  log  scales  where  the  floor  is  at  -60  dB.  This  is 
only  one-thousandth  the  height  of  the  peak  and  so  contributes  negligibly  to  the  image  MNR. 
However,  the  eye  views  image  intensity  logarithmically,  so  the  plot  was  scaled  as  such,  despite 
the  misleading  MNR  valu<;. 

The  spurious  targets  ate  also  present  in  the  ID  and  ID2  images,  and  there  is  little  that  can 
be  done  with  these  low-order  algorithms  to  reduce  their  height.  With  this  in  mind,  we  turn  to 


TABLE  5.2  Evaluation  Results  for  Fourier  Piece  (16,1)  on  Polar  Grid. 


Target 

Position 

Interp. 

Parameter 

MNR  (db) 

CPU  time 
(seconds) 

1  MNR/CPU  1 

-23.24 

NN 

-1.97 

3.28 

0.60 

G-ID 

(0.1) 

WBSM 

19.82 

Him 

G-ID 

(0.2) 

5.75 

ISIIM 

-23.24 

G-ID 

(0.3) 

mBBm 

22.57 

0.19 

-23.24 

wsinc 

2 

-4.13 

17.7 

miSMi 

-23,24 

wsinc 

4 

-7.45 

28.4 

-23.24 

wsinc 

6 

-10.62 

41.3 

0.26 

-23,24 

wsinc 

8 

-13.76 

53.0 

0.26 

-23.24 

wsinc 

10 

-17.01 

61.7 

0.28 

-23,24 

wsinc 

12 

-20.53 

70.6 

0.29 

-23.24 

wsinc 

14 

-24.48 

82.7 

0.30 

-23,24 

wsinc 

16 

-29.10 

94.7 

0.31 

-23.24 

wsinc 

18 

-34.66 

102.5 

0.34 

-23.24 

wsinc 

20 

-41.35 

112.0 

0.37 

-23.24 

wsinc 

22 

-47.10 

126.0 

0.37 

-23.24 

wsinc 

24 

-47.91 

134.0 

0.36 

-23.24 

B-spline 

-0.25 

-7.67 

pmmm 

-23.24 

B-spline 

-0.50 

-9.29 

Ktt 

-23.24 

B-spline 

-0.75 

-11.06 

17.22 

sKim 

-23,24 

B-spline 

-1.00 

17.22 

0.^6 

-23.24 

imn 

-15.35 

27.7 

0.55 

WlNOOW:  HAMHING  1 

MNR  I  IN  00)  :  -18.416^7 

FOURIER;  (  2.  8)  TARGET;  -23.30.  24.00  MAGITUOE;  1.0 


Figure  5.10  Inverse  Distance  Squared  Interpolator  at  (2.8). 
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the  more  sophisticated  weighted  sine  interpolators. 

It  was  shown  in  Chapter  3  that  the  sine  kernel  would  exactly  restore  a  bandlimited  signal 
which  was  properly  (Nyquist  rate,  infinite  length)  sampled.  The  weighted  sine  is  the  spatially 
limited  approximation  and  the  image  reconstructions  presented  here  demonstrate  that  even  a 
very  narrow  part  of  the  sine  performs  remarkably  well.  Figures  5.15  and  5.16  show  the  results 
of  interpolating  with  a  windowed  sine  of  orders  2  and  6  in  the  (2.8)  subarray.  Even  for  a  win¬ 
dow  of  only  2  samples  wide,  the  w-sinc  produces  an  MNR  of  -18.9  dB  which  is  much  better 
than  the  ID  interpolator.  The  CPU  time  is  even  lower  for  the  w-sinc  because  it  is  implemented 
separably.  If  the  window  is  widened  to  6  samples,  the  reconstruction  improves,  and  the  MNR 
drops  to  -30.4  dB.  The  w-sinc  size  can  be  increased  further  to  10  and  14  to  improve  the  MNR 
still  more  (Figs.  5.17  and  5.18),  but  soon  a  point  is  reached  where  the  MNR  fails  to  drop  (see 
Table  5.2).  The  reconstruction  is  approaching  the  exact  point  target  of  Fig.  5.2.  and  so  added 
terms  to  the  w-sinc  interpolator  simply  increase  the  algorithm  cost  without  a  noticeable  perfor¬ 
mance  improvement. 

If  the  w-sinc  interpolator  is  used  in  the  (16.1)  Fourier  subarray,  the  order  can  be  adjusted 
to  reduce  the  spurious  peak  to  an  acceptable  level.;  In  Fig.  5.19.  the  w-sinc  interpolator  of  order 
6  products  an  image  with  a  secondary  peak  at  -10.2  dB  below  the  peak.  If  the  order  is 
increased  to  16  (Fig.  5.20).  the  peak  drops  to  -28.6.  Of  course,  the  processing  time  rises  propor¬ 
tionately  from  41.3  seconds  to  94.7  seconds,  the  ratio  of  which,  0.44,  is  roughly  6/16  (0.38). 
i.e..  the  algorithm  processing  time  is  approximately  linear  with  interpolator  order.  If  the  order 
is  increased  to  20.  the  secondary  peak  disappears  from  the  image  (it  is  below  60  dB). 

The  processing  time  of  the  NN  and  GID  interpolators  is  relatively  small  compared  to  the 
w-sinc  of  any  order  greater  than  2.  The  strength  of  NN  and  GID  lies  in  speed,  though,  rather 
than  high  quality  image  reconstruction.  The  processing  times  for  each  of  the  various  interpola¬ 
tors  are  plotted  against  interpolator  order  in  Fig.  5.21.  The  NN  algorithm  is  referred  to  as 
zeroth  order,  and  the  GID  order  refers  to  the  power  of  the  inverse  distance.,  i.e..-  ID  is  first  order 
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Figure  5.19  W-Sinc  Interpolator  of  order  6  at  (16.1). 
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Figure  5.20  W-Sinc  Interpolator  of  order  16  at  ( 16,1 ). 
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and  ID2  is  second  order.. 

The  IMNRI/CPU  ratio  is  plotted  in  Fig.  5.22  lor  the  w-sinc  interpolator  as  applied  to 
several  of  the  subarrays.  Several  observations  can  be  made  by  studying  this  family  of  curves. 
First,  the  w-sinc  curves  all  converge  to  the  value  0.36  dli/s.  We  would  expect  convergence  to 
this  ratio  since  each  w-sinc  can  become  arbitrarily  close  to  the  exact  response  MNR.  and  the 
processing  time  of  w-sinc  is  independent  of  the  subarray.  Notice  that  the  ordering  of  the 
curves  from  top  to  bottom  corresponds  to  the  increasing  value  of  the  range  subscript  of  the 
subarray  coordinate  pair.  This  indicates  that  the  more  distant  subarrays  require  higher  order 
w-sinc  interpolators  to  achieve  the  same  MNR  values.  This  is  due  to  the  increasing  sample 
spacing  which  generates  the  spurious  targets.  These  targets  can  only  be  removed  by  increasing 
interpolator  order  at  the  expense  ol  processing  lime. 

The  IMNRI/CPU  value  can  be  thought  of  as  a  benefit/cost  ratio.  A  higher  ratio  indicates  a 
better  reconstruction  for  CPU  resources  used.  The  slope  indicates  the  amount  of  improved 
MNR  for  the  amount  of  additional  CPU  time  used.  Since  CPU  lime  is  proportional  to  interpola¬ 
tor  order,  a  more  negative  slope  corresponds  to  a  decrease  in  the  amount  of  MNR  improvement. 
Ultimately,  this  means  that  the  size  of  the  interpolator  can  be  varied  as  it  is  moved  out  in  the 
radial  direction  to  achieve  the  same  MNR  for  each  subarray.  This  concept  is  discussed  in 
Chapter  6  under  Further  Research. 

The  last  algorithm  shown  here  for  the  polar  grid  is  the  cubic  spline  and  B-spline  interpola¬ 
tors  described  in  the  last  section  of  Chapter  3.  The  first ' 'ersion  of  the  splines,  the  complete 
cubic  spline,  is  implemented  much  like  the  w-sinc  algorithm.  The  data  are  first  interpolated  to  a 
keystone  grid  by  calculating  the  spline  coefficients  for  each  interval  and  then  evaluating  the 
corresponding  cubic  polynomial  at  the  intermediate  (keystone)  points.  The  bulk  of  the  process¬ 
ing  is  consumed  in  calculating  the  coefficients  via  the  matrix  solution  step.  Evaluation  of  the 
cubic  polynomial  is  relatively  fast  compared  to  the  calculation  of  the  polynomial  coefficients. 
Next,  a  cubic  spline  is  generated  in  the  azimuth  direction  and  evaluated  along  the  rectangular 
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grid  coordinates.  As  shown  in  Table  5.1  for  the  subpatch  (2.8),.  the  MNR  is  roughly  that  of  a 
windowed  sine  of  an  order  between  6  and  8  (the  interpolator  order  need  not  be  even,  nor  must 
it  be  an  integer),  yet  processing  time  is  close  to  the  w-sinc  of  order  4.  The  cubic  spline  recon¬ 
struction  is  shown  in  Fig.  5.23.  The  IMNRI/CPU  ratio  is  much  higher  for  the  cubic  spline  than 
it  is  for  the  w-sinc.  since  it  performs  so  much  better  for  the  equivalent  processing  cost. 

The  modified  B-splines  also  do  quite  well  compared  to  the  w-sinc  interpolator.  They  were 
implemented  the  same  way  as  the  w-sinc.  i.e..  with  the  intermediate  keystone  grid.  Four 
different  values  of  the  parameter  in  Fq.  (3.52)  are  used.  -0.25.  -0.50,  -0.75.  and  -1.00.  with 
-1.0  providing  the  best  reconstruction.  The  MNR  of  -30.46  is  only  3  dB  worse  than  the  com¬ 
plete  cubic  spline,  and  the  processing  time  is  much  lower  (driving  the  IMNRI/CPU  ratio  up).  It 
is  shown  in  Fig.  5.24., 

Chapter  4  described  several  windows  which  may  be  applicable  for  SAR  data.  The  effects 
of  different  windows  may  be  seen  by  examining  the  image  reconstructions  from  a  typical  inter¬ 
polator:  for  example,  a  10th  order  w-sinc.  If  a  uniform  window  is  applied,  i.e..  no  weighting 
function,  then  the  reconstruction  is  a  very  narrow  spike  surrounded  by  some  low  sidelobes 
(Fig.  5.25a).  This  is  a  result  of  proper  target  placement  so  that  the  output  sine  is  sampled  in 
the  nulls.  The  sidelobes  are  a  result  of  interpolator  error.  If  the  Fourier  data  is  windowed  with 
a  disk  shaped  1-0  weighting  function,  then  the  MNR  is  dramatically  worsened  to  -4,87  dB  and 
the  image  looks  very  much  like  an  ,\'N  interpolation  (Fig.  5.25b). 

The  separable  Hamming  window,  the  standard  used  in  the  evaluations,  has  an  image 
reconstruction  of  the  same  data  shown  in  Fig.  5.26a.  The  MNR  is  a  very  low  -42.65  dB.  The 
output  of  a  circular  Hamming  window  is  much  noisier  in  appearance  (Fig.  26b)  and  has  an 
MNR  of  -28.12  dB.  It  is  apparent  that  the  discontinuity  at  the  Hamming  edge  is  causing  some 
problems.  There  is  an  additional.,  more  subtle  effect  occurring  as  a  result  of  using  the  circular 
Hamming  window.  The  zeros  of  the  circular  Hamming  Fourier  transform  do  not  fall  on  the 
Cartesian  grid  sample  points.  Thus,  even  an  exact  target  reconstruction  will  display  sidelobes 
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which  could  be  misleading.  This  can  be  corrected  slightly  through  the  use  of  the  Hanning  win¬ 
dow  which  goes  to  zero  smootnly.  The  Hanning  window  causes  faster  sidelobe  decay,  while 
having  a  slightly  higher  first  sidelobe  than  the  Hamming.  Figures  5.27a  and  5.27b  display  the 
results  of  the  Hanning  window,  both  separable  and  circular.  Note  that  the  separable  Hanning 
window  comes  very  close  to  matching  the  separable  Hamming  window  with  an  MNR  of  -40.59 
dB. 

5.5  Equi-PRF  Grid 

The  equi-PRF  grid  gained  little  in  the  way  of  performance..  Figures  5.28  and  5.29  are 
examples  of  the-  nearest  neighbor  algorithm  applied  to  the  equi-PRF  grid  in  the  regions  (2.8)  and 
(16,1)  which  may  be  compared  to  Figs.  5.5  and  5.8  for  the  standard  polar  format.  (Tables  C.5 
to  C.18  show  processing  time  and  MNR  values  for  each  of  the  equi-PRF  algorithms.)  The  pro¬ 
cessing  time  was  slightly  improved  for  the  w-sinc  algorithm  because  the  azimuth  pass  (interpo¬ 
lation  on  the  intermediate  keystone  grid)  calculations  were  simplified.  With  the  standard  polar 
grid,  the  keystone  samples  were  unevenly  spaced  in  azimuth,  resulting  in  excessive  coordinate 
positional  information  involving  trigonometric  functions.  With  equally  spaced  data,  positional 
information  is  calculated  as  a  simple  real  multiply..  For  splines  and  B-splines,  the  differences 
were  very  slight,  corresponding  to  the  time  taken  to  determine  into  which  interval  (between 
which  two  knots)  the  output  point  falls  prior  to  polynomial  evaluation.  With  non-uniform 
spacing  in  azimuth,  a  binary  search  is  performed  to  determine  the  interval,  but  with  equally 
spaced  data,  the  interval  is  again  found  with  one  real  multiply.; 

5.6  Keystone  Grid 

Since  the  keystone  grid  provides  data  along  vertical  lines,  the  interpolation  procedure  need 
act  only  in  one  dimension.  It  is  much  like  beginning  with  the  second  stage  of  the  w-sinc  out¬ 
put,  but  from  an  exact  first  stage  interpolation.  Thus,  a  better  reconstruction  should  be 
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SEPARABLE  RECTANGULAR  WINDOW 


m  :  -27. 43 


WINDOW:  RECTANGULAR  1 
ERIN  R5IZE:  10.0  ASIZE:  10.0 

FOURIER:  (  2.  8)  TARGET:  -23.00.  24.00  MAGITUDE:  1.0 
Figure  5.25a  Uniform  Window  After  10th  Order  W-Sinc  Interpolator. 
CIRCULAR  1-0  WINDOW 


nrm  :  -<1,87 
W1N00W:  RECTANGULAR  2 

ERIM  RS1ZE:  10.0  ASIZE:  10.0 

FOURIER:  (  2.  8)  TARGET:  -23.00.  24.00  MAGITUDE:  1.0 
Figure  5.25b  Circularly  Symmetric  Uniform  Window  (a  Disk). 
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POLAR:  WEIGHTED  SINC:  10 


HlNOOUi  HAMMING 
MNR  (IN  061  i  *42.64761 
INTERPs  ER1M  R-SlZEi 10.00  R-SIZE: 10. 00 
FOURIER:  t  2.  81  TARGE  it  -23.00.  24.80  MAGHUDE:  1.0 


Figure  5.26a  Interpolated  Data  With  Separable  Hamming  Window. 

’  CIRCULAR  HAMMING  WINDOW 


WlNOOW:  HAMMING  2 

ER1M  RS1ZE:  10.0  A5IZE:  10.0 

FOURIER:  '  (  2.  8)  TARGE1:  -23.00.  24.00  MAG  1 1 UOE :  1.0 


Figure  5.26b  Interpolated  Data  With  Circular  Hamming  Window. 
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EOUIPRFs  nearest  neighbor 


WNOOW:  HAMMING 

(IN  os,  .  .5.64513 

NEAREST  NEIGHBOR 

EGijR/ER:  (  2.  31  TARGET.  -23.00.  24.00  |  (J 

Figure  5.28  Equi-prf  Grid  wilh  Nearest  Neighbor  Interpolator  (2.8). 


E3UIPRF:  INVERSE  (0.3) 

I \ 


VNR  l IN  OBJ  :  -6.20923 
-  5IZE:  3  ORDER:  1,3 

COURIER:  (16.  91  TARGET:  -23.30.  2R.00  MAGNITUDE:  1.3 
Figure  5.29  Fqui-prf  Grid  with  Inverse  Distance  Interpolator  ( 16.1 ). 
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observed  (Figs.  5.30  through  5.41 ).  The  results  substantiate  this  theory  as  seen  in  Tables  C.9  to 
C.12.  livery  M\'R  value  is  noticeably  better  with  the  keystone  data  format  than  for  the  polar 
or  equi-PRF  formal.  It  is  interesting  to  note  that  the  w-sinc  interpolator  MNR  values  for  the 
keystone  grid  converge  much  faster  to  a  minimum  than  the  polar  format  data.  This  is  particu¬ 
larly  true  for  the  (16.1 )  subarray  where  the  spurious  targets  account  for  poor  \1\Rs.  Flimina- 
tion  of  the  range  line  interpolation  step  has  removed  the  error  introduced  into  the  azimuthal 
interpolation  stage,  in  the  polar  data  format,  this  error  could  only  be  minimized  by  increasing 
interpolator  order,  and  thus  processing  time,  to  obtain  the  same  MNR  as  the  keystone  data  for¬ 
mat.  The  interpolators  also  run  much  faster  because  the  spatial  position  calculations  are 
significantly  reduced  (only  one  dimensional  calculations)  and  the  two-stage  algorithms  now 
only  operate  in  one  stage. 

The  keystone  geometry  permits  the  use  of  the  chirp-z  algorithm  along  the  azimuth  data 
lines  as  described  in  Chapter  4.  This  combines  both  interpolator  and  FFT  sections  into  one 
stage.  Since  the  available  input  array  was  94  by  94  samples,  the  azimuth  input  vector  size  was 
94  with  an  output  vector  size  of  64  (in  the  spatial  domain).  The  results  of  the  chirp-z  algo¬ 
rithm  are  shown  in  Figs.  5.35  and  5.41  for  subarrays  (2.8)  and  (16.1).  respectively.  The  MNR 
values  of  -1.32  dB  for  (2.8)  is  surprisingly  poor  compared  to  the  other  keystone  interpolators, 
while  its  performance  in  the  (16,1)  region  was  similar  to  a  w-sinc  interpolator  with  order  12 
(in  speed  and  MNR),  The  processing  lime  is  high  because  the  FFT  included  in  the  algorithm 
must  be  zero  padded  to  perform  a  non-cyclic  convolution.  This  padding  increases  the  FFT  size 
from  64  to  256  (for  94  input  points.  256  is  the  first  power  of  two  greater  than  94+64), 

Note  in  Fig.  5.41  that  the  spurious  side  lobes  are  significantly  lower  than  the  w-sinc  order 
16  for  the  same  subarray.  This  is  because  the  chirp-z.  algorithm  interpolates/lransforms  in 
such  a  manner  that  the  spurious  target  analysis  of  Chapter  4  is  not  applicable.  The  high  spuri¬ 
ous  target  seen  in  Figs.  5.19  and  5.20  disappears  in  5.41,  because  the  sampling  rate  is  not 
changed  in  the  Fourier  domain,  but  rather,  the  transform  is  calculated  directly  from  the  input 
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SfUNE  I1HSL)  INOT-A-KNOT) 

FOURIER*  <  2.  8)  TARGET i  -23.18.  21.88  MAGNITUDE!  1.8 
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KEYSTONE:  WEIGHTED  SINC:  2 


MNR  UNOLI  i  -1.  36283 
HMGHTEO  SINC.  n-SlZEs  2.00 


FOURIERi  116.  II  IfWGtn  -23.00.  2^0*  "S»muK<  1.0 
Figure  5.38  Keystone  Grid  with  W-Sinc  Order  2  at  (16.1). 

KEYSTONE:  WEIGHTEO  SINC:  6 


Figure  5.39  Keystone  Grid  with  W-Sinc  Order  6  at  ( 16,1 ).. 
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data.  Theoretically.,  the  chirp-z  transform  should  behave  as  an  ideal  interpolator  (for  a  limited 
data  record)., 

5.7  Complexity  Analysis 

Using  the  methodology  for  complexity  analysis  discussed  above,  an  expression  for  each  of 
the  interpolators  was  calculated.  This  was  done  for  the  polar  and  keystone  grids.  Complexity 
measures  for  the  equi-prf  grid  are  similar  to  the  polar  format,  but  are  minus  the  trigonometry 
terms  which  are  required  for  the  data  position  calculations  in  the  second  stage  of  the  separable 
interpolators.  For  each  algorithm,  a  general  expression  is  presented  which  contains  a  break¬ 
down  of  the  software  in  terms  of  the  basic  cost  units  (Ca/s  .Cm/d,  Clrig,  etc.).  This  is  then 
reduced  to  an  intermediate  expression  by  making  certain  assumptions  about  the  input  array  size 
(number  of  input  azimuth  lines  and  N'2»N).  Finally,  a  simplified  expression  is  derived  by 
assuming  properties  of  the  data  processor  (computation  limes  for  a  C,rig.  Ca/s.  Cm/ri,  etc.).  It 
should  be  kept  in  mind  that  this  simplified  expression  is  only  as  accurate  as  the  assumptions 
made  in  the  approximation..  Recall  that  the  FFT  complexity  is  part  of  these  expressions.  It  can 
be  identified  as  the  2Nlog2N  term  in  the  expressions.  The  variables  N,  Ir.  Ia.  and  K  correspond 
to  the  output  grid  size  (N  by  N).  the  number  of  input  range  samples,  the  number  of  input 
azimuth  samples  and  the  interpolator  order  (w-sinc  only),  respectively.  Comparisons  of  the 
complexity  measures  are  done  with  N-1024  corresponding  to  the  original  image  size.  The 
nearest  neighbor  complexity  is  given  by  (5.1a)  through  (5.1c):, 

General  expression  for  nearest  neighbor  complexity: 

CNN  =  (5\2  +  N  +  2\2log2NX:iI1/(,  +  (5N2  +  N  +  2.\2log2.\)Cj/s  (5.1a) 

+  v?  ('  +  V*  c 

~  Vxtrig  *  ^  v'sqit 


Intermediate  expression  (N2»\): 
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CNN  =*  (5  +  2I.og2\)N2Clll/d  +  (5  +  21og2N)N2  Ca/S  (5.1b) 

+  \,2Cllls  +  N2Csl|n 

Simplified  expression  (Ca/s=Ctrif!=Cstir,=C.  Clll/d=5C): 

CNN  (32  +  121og2N)N2  C  (5.1c) 

Kquation  (5.1c)  shows  that  despite  the  simplicity  of  the  nearest  neighbor  algorithm,  it  is 
still  relatively  expensive  compared  to  the  FFT  times.  This  is  primarily  due  to  the  square  root 
and  trigonometric  calculations  required  to  locate  the  rectangular  output  points  in  polar  space. 

Since  the  ID2  interpolator  outperformed  ID  in  both  speed  and  MNK.,  it  was  the  only  ID 
type  of  algorithm  analyzed..  It  is  given  in  5.2a  through  5.2c., 

General  expression  for  Inverse  Distance  squared  (ID2)  complexi*... ' 

r1|ti  =  (32  +  2log2N)  N2  Cm/d  +  (29  +  2log2N)N2C#/s  +  N2Ctri)i 
Intermediate  expression  for  ID2  (la  =  lr  *  N) 

Crt)  **  (32  +  21og2N)N2Cm/d  +  (29  +  21og2N)N2Ca/s  +  N2Ctn| 

Simplified  expression  for  ID2  (Ca/s*C,rig*Csqr,=C.  Cm/d=5C) 

C,f,  (190+  121og2N)N2  C 

General  expression  Weighed  Sine  Interpolator  of  Order  K. 

Cwsmc  *  (41,  +  4.M,  +  3KNla  +  N  +  2Ma  +  3\2  +  2K  N2  +  2.\2log2N)Cm/d  (5.3a) 

+  (4la  +  5Nla  +  2kNIa  +  N  +  NIa  +  6N2  +  3KN2  +  2N2log2N)Ca/s 
+  (!a  +  Ma  +  2N,2)Ctl,g  +  (KNIa  +  KN2)Csinc 

Intermediate  expression  with  la  =  .N 

CWsmc=*(9  +  21og2N  +  5K)  N2  Cm/d  (5.3b) 

+  (12  +  21og2\  +  5K)\2  Ca/S 


(5.2a) 


(5.2b) 


(5.2c) 
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+  (3)  \-  Ctng 
+  2K  n2  rw 

Simplified  expression  (Ca/!>=CUl?=CStir,=C.  Cin/t,=5C) 

Cw«ntfe(62  +  12log2N  +  30k)  N2  C 


(5.3c) 


General  expression  for  Cubic  Spline  Interpolator: 

CS|>lme  =  (41,  +  10NIa  +  10NI,  +  N  +  8NI,  +  6N2  +  14  N2  +  2N2Iog2N)Cm/d  (5.4a) 
+  (4Ja  +  10NI,  +  16N1,  +  2N  +  4Nla  +  ION2  +  20N2  +  2N2log2N)Ca,s+  (Ia  +  N)  CUlR 


Intermediate  expression: 


+  2logA)  N-  CmM 

+  l60+21og2N)N2C„, 

+  2N  CM1K 


(5.4b) 


Simplified  expression: 


CS|»lmeSs(300  +  121og2N)  N2  C 


(5.4c) 


The  cubic  spline  interpolator  has  a  very  high  order  of  complexity  due  to  the  matrix  solution 
step.  It  has  an  interpolator  complexity  which  approximates  that  of  the  w-sinc  of  order  S. 
Tables  C.l  through  C.4  show  that  the  performance  is  similar  to  the  w-sinc  with  an  order 
between  6  and  8. 


General  expression  for  B-spline  Interpolator:. 

CB— s,>ii«ie  =  (41,  +  4Ma  +  12M,  +  N  +  2NIa  +  3N2  +  H  N2  +  2\2Iog,\)Cm/<i  (5.5a) 

+  (41,  +  5M,  +  8Ma  +  N  +  Ma  +  6N2  +  12N2  +  2N2log2N)Ca/a+  (Ia  +  Nla)  Ctllt, 

Intermediate  expression: 


CB_sl)hlK==(29  +  21og3\)N2Cl 


iii/d 


(5.5b) 
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+  (30  +  21og2N)  N2  Ca/s 

+  i\2C,flg 

Simplified  expression:. 

Cb-s,,,.,*5^  176  +  1 21og2N)  N2  C  (5.5c) 

The  13-splines  have  a  complexity  which  is  approximately  half  that  of  the  complete  cubic 
spline,  yet  their  performance  is  comparable.  The  complexity  also  approximates  the  w-sinc  of 
order  4.  This  is  expected,  since  the  fourth  order  w-sinc  and  B-spline  interpolators  operate  the 
same  way:  by  convolving  with  4  input  points  in  2  stages. 

The  keystone  grid  reduced  the  complexity  measure  by  roughly  a  factor  of  2.  This  is  due 
to  the  elimination  of  one  stage  of  the  separable  interpolators  and  the  2D  positional  calculations 
for  the  lower  order  algorithms.  Equations  (5.6)  through  (5.8)  present  the  complexity  expres¬ 
sions  for  the  interpolators  used  on  the  keystone  grid.. 

General  expression  for  complexity  of  NN  on  keystone  grid: 

«  (3.\2  +  2N2log2N  +  N)Cm/<i  +  (3N2  +  2N2log2N  +  N)Ca/s  +  n2C,rig  (5.6a) 
Intermediate  expression: 

Cfff  **  (3  +  21og2N)N2Cm/d  +  (3  +  21og,N)N2Ca/s  +  N2Ctrig  (5.6b) 

Simplified  expression: 

C$f  a®  ( 1 9  +  1 2  log2N )  N2C  (5.6c) 

General  expression  for  linear  interpolator  on  keystone  grid: 

C,*'>ar  =  (Ia  +  31a!r  +  N  +  15N2  +  2N2log2N)Cl,1/d  (5.7a) 

+  (Ia  +  lalr  +  N  +  14N2  +  2N2log,N)Ca/s 
+  (21a  +  N2)Ctrig 


Intermediate  expression: 
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Cfer  =  (18  +  21og2N)N2Cm/d  +  (15  +  21og2N)N2Ca/s  +  N2Ctnf, 
Simplified  expression:' 


Qmenr  28  (1«»6  +  l21og2.\)N2C 


(5.7b) 


(5.7c) 


General  expression  for  weighted  sine  on  keystone  grid:; 

C*gie  =  (N  +  2Nla  +  3!\2  +  2N2K  +  2N2log2N)Cm/(1  (5.8a) 

+  (N  +  Ma  +  6N2  +  3N2K  +  2N2log2N)Ca/s 
+  (Nla  +  2Na)Ctt„  +  N2KCSII)C 

Intermediate  expression:. 

C&L  =b  (5  +  2log;N  +  2K )N2Cm/,i  +  (7  +  21og2N  +  2K)N2Ca/s  (5.8b) 

+  3\2C(ne  +  N2KCimt 

Simplified  expression:' 

C&L  28  (35  +  121og2N  +  13K)N2C  (5.8c) 

Using  Eqs.  (5.6c).  (5.7c).  and  (5.8c)  with  N»1024.  and  K-12.  it  is  found  that 
28  0.91  CNN.  C&  0.72  C,i.  and  C*I'smc  0.57  Cw_imc.  The  complexity  of  the  NN  on 
the  keystone  grid  has  not  improved  much  over  the  polar  grid.  This  is  because  the  FFT  com¬ 
plexity  is  dominating  both  expressions.  The  \x-sinc  interpolator,  which  is  dominated  by  the 
interpolation  step,  shows  a  noticeable  cost  improvement,  owing  to  the  removal  of  the  range  line 
interpolations. 


General  expression  for  the  cubic  spline  on  keystone  grid: 

C*|  L  =  (N  +  3N2  +  8X2  +  2N2log,N)CmA,  (5.9a) 

+  (6.N  +  27M,  +  8N2  +  2\2log2.\X-a/s 
+  I  C 

~  »a  Mug 

Intermediate  expression: 

C4.1L  =  (11  +  21og2\).\2Clll/d  +  (35  +  21og2\  )N2CJ/S 


(5.9b) 
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5.8  Summary 

This  chapter  has  presen  let!  several  evaluations  of  the  interpolators  described  and  analyzed 
in  Chapter  4.  The  required  interpolator  accuracy  is  dependent  upon  the  acceptable  reconstruc¬ 
tion  accuracy,  measured  here  as  a  mulliplicative-noise-ralio,  and  on  the  subarray  location 
within  the  toroidal  slice.  The  more  radially  distant  subarrays  require  a  higher  order  algorithm 
to  reduce  the  MNR,  sidelobes,  and  spurious  targets,  while  the  subarrays  closer  to  the  origin  are 
nearly  rectangular  and  require  relatively  little  computation  to  produce  a  good  reconstruction. 
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CHAPTER  6 

CONCLUSIONS  AND  FURTHER  RESEARCH 

The  problem  of  SAR  data  interpolation  in  the  Fourier  domain  was  examined,  and  several 
types  of  interpolator  algorithms  were  discussed,  analyzed,  and  evaluated  experimentally.  A 
true  analysis  of  the  polar  to  rectangular  interpolation  problem  seems  intractable,  so  the  prob¬ 
lem  was  reduced  to  a  rectangular-rectangular  problem  and  finally  a  one-dimensional  problem 
of  interpolation  in  the  Fourier  domain.  Originally,  much  work  was  done  to  simulate  an  exact 
target  response,  but  this  obscured  the  effects  that  were  being  studied.  The  entire  SAR  and  tomo¬ 
graphic  grid  interpolation  problem  is  that  of  finding  an  interpolator  that  is  reasonably  fast  and 
vet  has  a  transform  that  is  close  to  an  ideal  low-pass  filler.  The  various  two-dimensional  inter¬ 
polators  that  were  in  current  use  [2\]  were  more  carefully  studied,  and  some  new  interpolators 
for  Fourier  space  were  proposed  for  the  SAR  problem,  i.e.,  inverse  distance  squared,  windowed 
sine,  and  cubic  splines. 

The  inverse  distance  squared  interpolator  has  not  been  seen  in  the  recent  DSP  literature, 
but  proved  to  generate  good  reconstruction  at  a  lower  cost  than  the  more  heavily  used  inverse 
distance  algorithm.  All  of  the  known  examples  which  examine  DSP  interpolators  compare  a  new 
kernel  with  a  sine  or  truncated  sine.  It  seems  that  this  comparison  is  unfair,  since  the  win¬ 
dowed  sine  has  a  much  better  response  and  is  not  too  difficult  to  compute.  Spline  interpolation 
has  also  been  successfully  used  in  the  area  of  Fourier  domain  reconstruction.  This  work  helps 
sort  out  the  meaning  of  spline  in  the  current  literature. 

It  is  unlikely  that  an  interpolation  kernel  will  be  discovered  which  is  both  easy  to  imple¬ 
ment  (fast)  and  has  the  desired  spectrum  (ideal  low-pass  filter)  since  the  two  criteria  work 
against  each  other.  Classical  optimization  procedures  are  ineffective  here  because  the  algorithm 
cost  function  is  too  difficult  to  paramelize  well.  The  weighted  sine  is  the  best  example  of  an 
interpolation  algorithm  whose  order  can  be  adjusted  to  reduce  interpolation  error  below  a  given 
specification,  though  at  the  expense  of  processing  time.  The  optimal  Fourier  domain  interpolator 


140 


is  global  in  nature,  taking  all  the  known  data  into  account.  However,  this  usually  leads  to 
prohibitive  memory  usage  and  computation  times.  Better  methods  exist  for  producing  better 
image  quality,  such  as  iterative  techniques,  but  these  are  far  and  away  too  expensive  for  the 
type  of  real-time  processing  that  is  desired. 

A  new  method  of  evaluating  the  windowed  point  target  response  was  presented  as  an 
alternative  to  the  MSI:  criterion  often  used  in  image  processing  -  the  MNR  and  MNR/CPU  ratio. 
These  figures  of  merit  proved  useful  in  rating  the  interpolator  image  quality  and  in  making 
algorithm  comparisons. 

A  novel  approach  to  the  interpolation-inversion  stage  was  presented  via  the  chirp-z 
transform.  Although  it  did  not  seem  competitive  with  the  weighted  sine  interpolator,  it  has 
promise  as  a  good  alternative.. 

The  alternative  sampling  grids  reduced  the  interpolation  error  dramatically  and  were 
much  faster  to  implement,  due  to  the  one-dimensional  nature  of  the  reconstruction.  It  is  sug¬ 
gested  that  these  raster  designs  be  used  in  actual  hardware  designs  due  to  the  tremendous  CPU 
cost  savings  they  offer. 

In  short,  the  algorithm  which  produced  the  worst  reconstruction  was  the  nearest  neighbor,; 
followed  by  the  inverse  distance  squared,  and  then  inverse  distance.  These  algorithms  were  not 
competitive  with  the  separable  interpolators:  the  weighted  sine  and  cubic  splines.  The  weighted 
sine  had  the  advantage  of  having  an  adjustable  parameter  (order)  which  could  be  increased  to 
improve  reconstruction  to  acceptable  levels,  while  the  cubic  spline  was  easier  to  implement  for 
the  same  level  of  image  quality.  The  chirp  z-transform  produced  a  good  reconstructed  image, 
yet  was  not  very  cost  competitive  with  the  separable  algorithms  due  to  the  extended  FIT  size 
in  the  convolution  stage. 
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6.1  Further  Work 

The  advent  of  digital  processing  of  SAR  signals  has  opened  a  whole  area  of  research.  The 
interpolation  problem  examined  in  this  work  is  useful  in  the  direct  inversion  of  the  Fourier 
data  set:  however,  there  are  other  inversion  algorithms  which  may  lead  to  faster  reconstruc¬ 
tions  if  they  in  turn  are  made  faster,  specifically,  convolutional  back-projection.  The  work  here 
has  also  led  to  many  other  questions  concerning  2D  interpolation  and  SAR.  Additional  research 
topics  are  presented  below, 

6.1.1  Full  data  array  evaluations 

The  computer  evaluations  in  this  work  were  done  on  a  relatively  small  data  set.  This  had 
the  advantage  that  the  effects  of  the  interpolator  in  different  Fourier  regions  could  be  studied. 
However,  an  actual  implementation  of  a  1024  by  1024  (or  larger)  data  array  would  generate  an 
image  which  is  the  coherent  sum  of  all  the  small  sub-arrays.  It  would  be  very  useful  to  work 
with  the  large  polar  grid  interpolation  problem  and  examine  the  results  for  any  new 
phenomena  which  may  appear.  The  use  of  the  newer  generation  supercomputers  seems  ideal 
for  this  higher  order  problem, 

6.1.2  Oversampling 

As  mentioned  in  Chapter  4,  the  true  input  data  of  a  working  system  is  sampled  at  a  rate 
much  higher  than  required  for  the  given  system  resolution.  Work  should  be  done  to  see  how 
this  higher  volume  of  data,  and  hence  finer  sample  spacing,  could  be  used  prior  to  prefiltering  to 
reduce  interpolation  error. 

6.1.3  Spline  approximation  to  the  sine 

Splines  have  been  increasingly  popular  in  DSP.  Though  the  usual  Fourier  domain  interpre¬ 
tation  is  lacking,  they  are  easy  to  evaluate  and  have  some  nice  mathematical  properties.  Since  it 
was  shown  that  the  weighted  sine  provided  the  best  reconstruction,  it  may  be  useful  to  see  if  a 
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cubic  spline  approximation  to  the  sine  would  prove  to  be  as  accurate. 

6.1.4  Spatially  varying  interpolation  order 

It  was  seen  that  different  parts  of  the  Fourier  data  set  require  different  degrees  of  interpo¬ 
lation  accuracy  based  on  the  sample  spacing  and  angular  orientation  to  the  output  grid.  It  is 
suggested  that  the  interpolator  order  be  made  spatially  dependent  to  minimize  the  amount  of 
computation  required,  i.e.,  a  low-order  sine  could  be  used  in  the  (1.1)  region  of  the  torus  and  a 
higher  order  interpolator  could  be  used  in  the  more  (radially)  distant  parts  of  the  polar  grid 
(16.fi).  Kmpirically.  it  seems  that  the  required  interpolator  order  is  more  dependent  on  sample 
rate  changes,  resulting  in  spurious  targets,  than  the  degree  of  grid  rotation,. 

The  order  of  the  separable  sine  interpolator  studied  here  was  always  the  same  in  both 
radial  and  azimuthal  directions.  This  may  not  be  needed,  as  the  azimuthal  rate  change  demands 
a  higher  order  interpolator  than  the  range  lines,  further  computational  savings  may  result  if 
the  range  and  azimuth  interpolator  orders  are  minimized  to  correspond  to  a  prescribed  accuracy 
for  each  dimension. 
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APPENDIX  A 

A  FAST  EVALUATION  OF  A  PERIODICALLY  SAMPLED  SINE 

A  major  obstacle  in  the  use  of  the  weighted  sine  kernel  is  the  cost  of  evaluating  the  two 
transcendental  within  the  kernel  function: 

(.54+.46cos(x*c2)  )*sin(x*Cj) 

hU)  =  - ; -  (A.l) 

x*ct 

where  x/cj  =  rr  and  c2  determines  the  extent  of  the  Hamming  window  (the  interpolator  order). 

This  can  be  approached  in  many  ways.  First  is  the  direct  evaluation  method.  The  inter¬ 
polation  kernel  g(x)  is  a  simple  evaluation  using  floating  point  sine  and  cosine  libraries.  While 
this  produces  very  accurate  results,  it  is  the  most  time-consuming,  because  the  transcendental 
are  usually  computed  with  a  power  series  expansion  with  many  terms. 

An  alternative  to  the  direct  method  is  with  a  table  lookup.  Here,  the  sine  function  is 
stored  as  a  finely  sampled  array  stored  in  a  ROM.  and  the  values  of  the  transcendental  are  cal¬ 
culated  by  finding  the  closest  value  in  the  table.  A  finer  evaluation  may  be  obtained  by  linearly 
interpolating  between  the  two  nearest  table  entries.  This,  however,  can  produce  additional  error 
in  the  kernel.  Also,  the  error  in  the  sine  function  greatly  increases  as  x  approaches  zero.  This  is. 
of  course,  due  to  sin(x)  and  x  approaching  zero  at  the  same  rate.  Machine  precision  begins  to 
cause  errors  in  the  quotient.  To  correct  this,  we  can  store  the  sinc(x)  as  one  table  and  the 
weighted  cos(x)  as  another.  Note.  loo.  that  for  storage  savings,  only  the  positive  half  of  sinc(x) 
needs  to  be  stored,  and  only  one-fourth  of  the  cosine  function  must  be  tabled.  Half-angle  for¬ 
mulas  may  be  used  to  index  to  the  correct  table  value. 

If  the  data  are  equally  spaced,  a  third,  recursive,  method  may  be  used.  Because  we  are 
computing  a  sinusoid  at  constantly  spaced  points,  we  can  take  advantage  of  the  complex 
exponential  to  recursively  yield  successive  sine  values.  We  wish  to  calculate  sin(n*0  +  a)  with 
n  being  the  order  of  the  interpolator.  6  and  a  are  determined  by  the  position  of  the  kernel  func¬ 
tion  relative  to  the  known  data.  We  know  that 


r- 
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ImagleJ0)  =  sin(a) 

Using  this  identity,  we  can  compute 

(A. 2) 

S„  =  sin(n*0  +  a)  For  n  =  0 . N— 1 

(A. 3) 

sin(n*0  +  a)  =  lmagleJ(,>-f>  +  o)} 

(A. 4a) 

=  Imagle^V*} 

Let 

(A.4b) 

«nrv 

II 

(A.5) 

and 

or,,  =  e^ 

(A.6) 

Then  or,  is  recursively  defined  as 

o-.-oVi  *£  i  =  1....N-1 

Then: 

(A. 7) 

S,  *  Imagid,}  i  *  1 . N— 1 

(A. 8) 

This  is  a  significant  computational  improvement  over  the  standard  power  series  expansion 
which  is  computed  for  each  evaluation.  The  cosine  term  of  A. I  may  be  computed  in  the  same 
manner. 


\ 
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APPENDIX  B 

PLACING  THE  DC  POINT  IN  THE  IMAGE  CENTER 

The  final  stage  of  the  SAR  processing  includes  a  circular  permutation  of  the  2D  data  set  to 
place  the  dc  point  in  the  center..  Since  this  is  actually  the  spatial  domain,  this  corresponds  to 
the  terrain  center.  This  step  may  also  be  performed  in  the  Fourier  domain  with  a  slight  compu¬ 
tational  savings. 

As  shown  by  Rosenfeld  and  Kak  [74]  the  output  of  the  FFT  will  be  rotated  by  one-half 
the  image  size  in  both  coordinates  if  it  is  first  multiplied  in  the  frequency  domain  by  (— l)m+n. 
That  is.  given  the  original  Fourier  data  array  F(m.n).  define  F(m.n): 

F(m.n)  *  F(m.nM— l)m+n  (B.l) 

A 

then  the  inverse  transform  of  F(m.n)  becomes 


-  M-«  J2*  ^  ♦  £? 

f(j.k)  *  £  £F(m.nM-l)m+ne  M 

(B.2) 

j*o  k«0 

but  we  can  rewrite  (— l)nH'n  as 

_  ,  m  +  n 

(-l)ro-Hi  «  e  **  ~TT 

(B.3a) 

m  £-pw  (m/2  ♦  n/2) 

(B.3b) 

-  (M/2)m  (N/2  )n 

M  HKT~ 

*  e 

(B.3c) 

substituting  (B.3c)  into  (B.2)  and  gathering  terms  yields 

f(j.k)  =  £  £F(m.n)e 

(B.4) 

-  f(  j — M/2,  k— N/2)  (QED) 

T 
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APPENDIX  C 

EVALUATION  RESULTS  FOR  THREE  GRID  FORMATS 

TABLE  C.l.  Evaluation  Results  for  Fourier  Piece  (2.8)  on  Polar  Grid, 


Target 

Position 

Interp. 

Parameter 

MNR  (dB) 

CPU  time 
(seconds) 

I  MNR/CPU  1 

-23.24 

NN 

-5.64 

3.28 

1.72 

-23.24 

G-ID 

(0.1) 

-14.56 

19.82 

0.73 

-23.24 

G-ID 

(0.2) 

-18.41 

5.75 

3.20 

-23.24 

G-ID 

(0.3) 

-17.00 

22.57 

0.75 

-23.24 

G-ID 

(0.4) 

-14.94 

21.17 

0.71 

-23.24 

G-ID 

(0.5) 

-13.37 

19.37 

0.69 

-23.24 

wsinc 

2 

-18.92 

17.7 

1.07 

-23.24 

wsinc 

4 

-24.50 

28.4 

0.86 

-23.24 

wsinc 

6 

-30.37 

41.3 

0.74 

-23.24 

wsinc 

8 

-36.53 

53.0 

0.69 

-23.24 

wsinc 

10 

-42.65 

61.7 

0.69 

-23.24 

wsinc 

12 

-46.88 

70.6 

0.66 

-23.24 

wsinc 

14 

-48.06 

82.7 

0.47 

-23.24 

wsinc 

16 

-48.24 

94.7 

0.51 

-23.24 

wsinc 

18 

-48.27 

102.5 

0.47 

-23.24 

wsinc 

20 

-48.27 

112.9 

0.43 

-23.24 

wsinc 

22 

-48.30 

126.0 

0.38 

-23.24 

wsinc 

24 

-48.38 

134.0 

0.36 

-23.24 

B-spline 

-0.25 

-23.01 

17.22 

1.34 

-23.24 

B-spline 

-0.50 

-25.17 

17.22 

1.46 

-23.24 

B-spline 

-0.75 

-27.62 

17.22 

1.60 

-23.24 

B-spline 

-1.00 

-30.46 

17.22 

1.77 

-23.24 

Spline  (IMSL) 

-32.41 

27.7 

1.17 

TABLE  C.2.  Evaluation  Results  for  Fourier  Piece  (10.5)  on  Polar  Grid. 


TABLE  C.3.  Evaluation  Results  for  Fourier  Piece  (16.1)  on  Polar  Grid. 


TABLE  C.4.  Evaluation  Results  for  Fourier  Piece  (16,8)  on  Polar  Grid. 


TABLE  C.5.  Evaluation  Results  for  Fourier  Piece  (2.8)  on  Equi-PRF  grid. 


TABLE  C.7.  Evaluation  Results  for  Fourier  Piece  (16.1)  on  Equi-PRF  grid. 
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TABLE  C.8.  Evaluation  Results  for  Fourier  Piece  (16.8)  on  Equi-PRF  grid. 


Target 

Position 

Interp. 

Parameter 

MNR  (dB) 

CPU  time 
(seconds) 

1 MNR/CPU  1 

-23,24 

NN 

-1.00 

3.58 

0.28 

-23.24 

G-ID 

(0.1) 

-6.21 

19.02 

0.326 

-23.24 

G-ID 

(0.2) 

-4.97 

5.67 

0.88 

-23.24 

G-ID 

(0.3) 

-4.55 

18.75 

0.24 

-23.24 

G-ID 

(1.1) 

-2.98 

63.62 

0.047 

-23.24 

'  G-ID 

(1.2) 

-3.38 

17.15 

0.20 

-23.24 

G-ID 

(1.3) 

-3.90 

70.37 

0.055 

-23.24 

G-ID 

(2.1) 

-1.28 

151.82 

0.008 

-23.24 

G-ID 

(2.2) 

-2.95 

20.48 

0.144 

-23.24 

G-ID 

(2.3) 

-3.79 

207.55 

0.018 

-23.24 

wsinc 

2 

-4.04 

11.98 

0.337 

-23.24 

wsinc 

4 

-7.62 

22.57 

0.338 

-23.24 

wsinc 

6 

-10.81 

32.70 

0.331 

-23.24 

wsinc 

8 

-13.89 

41.10 

0.338 

-23.24 

wsinc 

10 

-17.07 

50.73 

0.336 

-23.24 

wsinc 

12 

-20.55 

59.78 

0.344 

-23.24 

wsinc 

14 

-24.52 

69.53 

0.353 

-23.24 

wsinc 

16 

-29.24 

79.20 

0.369 

-23.24 

wsinc 

18 

-35.04 

88.93 

0.394 

-23.24 

wsinc 

20 

-42.30 

99.03 

0.427 

-23.24 

wsinc 

22 

-48.06 

1 10.03 

0.437 

-23.24 

wsinc 

24 

-47.44 

122.02 

0.389 

-23.24 

B-spline 

-0.25 

-8.32 

22.13 

0.38 

-23.24 

B-spline 

-0.50 

-9.03 

22.13 

0.41 

-23.24 

B-spline 

-0.75 

-11.83 

22.13 

0.54 

-23.24 

B-spline 

-1.00 

-14.33 

22.13 

0.65 

-23.24 

Spline  (IMSL) 

-42.57 

26.75 

1.59 

TABLE  C.9.  Evaluation  Results  for  Fourier  Piece  (2.8)  on  Keystone  Grid. 


Target 

Position 

Interp. 

Parameter 

MNR  (dB) 

CPU  time 
(seconds) 

i  MNR/CPU  1 

-23.24 


-23.24 

-23.24 


23.2 

23.2 


-23.24 

-23.24 

-23.24 

-23.24 

-23.24 

-23.24 

-23.24 


-23.24 


23.2 

23.2 


-23.24 

-23.24 


-23.24 


NN 


G-1D 

G-1D 

G-ID 

G-1D 


wsinc 

wsinc 

wsinc 

wsinc 

wsinc 

wsinc 

wsinc 

wsinc 

wsinc 


Chirp-Z 


B-spline 


Spline  (IMSL) 


-7.26 


-24.16 

-22.23 

-17.66 

-15.15 


-22.60 

-26.95 

-31.86 

-37.43 

-43.14 

-46.97 

-48.04 

-48.12 

-48.14 


-1.32 


-34.12 

-35.88 

-36.12 

-36.98 


-34.80 


0.67 


0.85 

1.42 

1.65 

1.47 


4.63 

8.70 

13.2 

17.3 
21.90 
26.00 
29.85 
34.52 
37.57 


26.50 


.95 
.95 
7.95 
7.95 


12.06 


10.83 


10.70 

10.31 


4.88 

3.10 

2.41 

2.16 

1.97 

1.81 

1.61 

1.39 

1.28 


0.05 


.29 
.51 
.54 
4.65 


2.89 


TABLE  C.10.  Evaluation  Results  for  Fourier  Piece  (16,1)  on  Keystone  Grid. 


TABLE  C.ll.  Evaluation  Results  for  Fourier  Piece  (10,5)  on  Keystone  Grid. 
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